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ABSTRACT 


Travel time of an acoustic signal from transmitter to receiver provides a 
great deal of information about the ocean environment. Variations in the trav¢’ 
time of the signal may be caused by the changes in the sound speed along the 
path. Since sound: speed is a function of pressure, temperature and salinity, 
measurement of this parameter in acoustic tomography provides a means to 
ebserve ocean fluctuations through the use of inverse techniques. The upcoming 
Heard Island Experiment will attempt to determine the feasibility of measuring 
global warming by measuring changes in signal travel time that may be caused 
by temperature changes in the world's oceans. The signals to be transmitted in 
this experiment are phase encoded maximal-length sequences of various lengths 
which are well suited to measurement of travel time. The objectives of this 
thesis are to provide a software package, in C, that will allow participation as a 
receiver in this experiment, and to provide a general. capability to process any 
maximal-length sequence, transmitted at any carrier frequency and with any 
reasonable Doppler. A background on wave propagation, maximal-length 


sequences, and Doppler processing are presented in this thesis. 
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I. INTRODUCTION 


A. THESIS SUMMARY 

The.objective of this thesis is to develop a program in the C programming 
language [Ref. 1] that can process acoustic signals used in ocean acoustic 
tomography. The goal of tomography signal processing is the precise 
measurement of acoustic travel time, the integral along the raypath of inverse 
sound speed. Travel time is a fairly well understood function of temperature, 
salinity, and pressure. Once variations in the travel time are known, as well as 
the travel times for multiple arrival. paths, the fluctuations of the ocean can be 
determined from mathematical inverse techniques [Ref. 2]. 

One method in which travel time can be measured is through the use of 
explosive and implosive devices. These crude tools, however, are not exactly 
repeatable. A better method that has been found employs the use of maximal- 
length sequences as a phase modulation for bandpass signals. Maximal-length 
sequences (m-sequences) or pseudo-random noise are well suited for this 
application because of their deterministic nature, correlation properties, and 
simplicity. 

The goal of the work described in this thesis was to be able to measure 
arrival time by performing replica cross-correlation of a transmitted acoustic 
signal that has been p! ase-encoded using m-sequences, and to be able to detect 


that signal over any reasonable Doppler shift. Specifically, the programming 


package should be able to: 




















1. Detect any phase encoded m-sequence transmitted with any carrier 


frequency. 

ae Scan over any reasonable Doppler range at any Doppler bin 
spacing. 

3. Provide filtering capability with any type filter. 

4. Employ fast algorithms to minimize processing time. 

Se Be well documented and portable for future transfer to a dedicated 


machine for real time processing. 
The body of this thesis is structured as follows: 
Chapter IT]. Acoustic Wave Propagation. 
Chapter II. M-Sequences and Hadamard Processing. 
Chapter IV. Time-Varying Doppler -Processing. 
Chapter V. Results and Conclusions. 

Chapter {fT presents an introduction to acoustic wave propagation in the ocean 
environment. It includes plane wave propagation, spherical wave propagation, 
the dependencies of sound speed on temperature and other parameters. Acoustic 
raypath determination is discussed. 

The third chapter is an introduction to general shift registers and m- 
sequences. The Fast Hadamard Transform is reviewed to provide the 
background necessary to understand the programming package. 

Chapter IV presents the signal processing aspects of the thesis. It first 
presents basic phese encoding using m-sequences, and then develops the treatment 
of a signal with zero Doppler. And finally, the nonzero Doppler case is 
presented. 

Appendix A contains the results from all Doppler considerations. Appendix 


B is the source code for the thesis work as well as some supplementary 











programs. Supplementary programs include: plotting routines using VAX NCAR 
Graphics {Ref. 3], a routine for generating shift register states and corresponding 
m-sequences in a (1,-1) format for viewing before processing, if desired. 
Finally, a MATLAB program that can be used for generating filter coefficients 
for a Butterworth bandpass filter and a Chebychev lowpass filter. The filter 
coefficients are stored to a file that can be read by the main program SEQREM. 
Documentation of the sorrce code should be sufficient for the typical-user to 
understand the general operation of the program without much difficulty, and is 
built using standard C functions and coding to minimize any portability 


problems. 


B. THE HEARD ISLAND EXPERIMENT 

The first direct application of this thesis is planned for the beginning of 
1991. From January 23 - February 4, 1991 an experiment to determine the 
feasibility of measuring global warming using underwater acoustic signals is to 
take place and is currently designated the Heard Island Experiment (Ref. 4]. It is 
hoped that changes in the temperature of the Earth's ecosystem can be measured 
by observing changes in travel time of an acoustic signal over distances which 
include one or more oceans. Travel time is related to water temperature, as well 
as other parameters, so that a change in travel time might be related to warming 
of the Earth's atmosphere, which would also cause a overall warming of the 
world's oceans. If these two parameters can be related, then existence of global 
warming might be determined. 

The signals will be transmitted from the vicinity of Heard Island in the 
southern Indian Ocean with raypaths that reach several receiving sites in the 


Indian Ocean, the Adantic, and the Pacific as shown in Figure 1.]. The 











transmitting vessel will be the m/v Cory Chouest. Preliminary ray traces have 
shown that one possible reception site is the Monterey Peninsula [Ref. 5]. 
Participation is intended by the Naval Postgraduate School, and is one motivation 
for this thesis. Of the many signals that will be transmitted, one is an m- 
sequence of 255 digits with a Q of 5, and a carrier frequency of 57 Hz, where Q 
is the ratio of digit frequency to carrier frequency and is normally selected to be 
an integer number. A carrier of 57 Hz was chosen because of its long 
propagation distance and also to help distinguish it from the common frequencies 
50 and 60 Hz that are generated by a large number of power plants and 
machinery world wide. 

Detailed signal specifications are contained in Reference 4, but some of the 


more general data are: 


1. HLF4LL source. 
2s Output power: 209 dB. 
3. Signal to be transmitted using m-sequences 


s(t) = Acos(2nft + M(t)y) 


where f, = 57 Hz. 





M(t) is the m-sequence. 
w= 45 degrees, is the modulation angle. 


4, Maximum of 3-5 knots to maintain way. 
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II. ACOUSTIC WAVE PROPAGATION 


A. THE WAVE EQUATIONS 

Ocean acoustic signals frequently travel from source to receiver along 
multiple paths. These multiple "raypaths" are formed as a direct result of the 
nonhomogeneous structure of the ocean. However, understanding how waves 
propagate in the ocean can be facilitated by first treating it as a homogeneous 
medium and developing equations for wave propagation, and next, by the 
introduction of Snell's Law and some basic rules for determining the raypaths in 
a nonhomogeneous medium. 

Sound travels as a result of the displacement of particles at some source, 
which causes a local change in pressure. The resulting pressure change then 
propagates away from the source in a spherical manner as a pressure wave with 
velocity c. At distances sufficiently far from the source the wave can be treated 
as planar over a finite region. 

1. Plane Wave Equation 

Assuming the propagation direction to be in the x direction with the 
wave front formed in the y-z plane, the pressure change and particle speed can 
be related by the equation [Refs. 6,7,8] 


dp(x,t) % pole 


2.1 
ox or’ oo 


where p(x,t) is pressure, u(x,t) is particle speed, and p is the density of water. 


The rate of change of u(x,t) with respect to x is given by [Ref. 6] 








duis) _ 1 apost) 


22 
ox B at’ Co 
. where B is the bulk modulus of elasticity. Combining Eq. (2.1) and Eq. (2.2) 
gives 
a, (x,t) _ B 2 (x,t) 
p(x,t p(x,t 
mE og. Teh. Faye a (2,3) 
p 
at Ox 
The sound speed in water, c, is related to B and p and is given by 
a3, (2.4) 


p 
Substituting into Eq. (2.3) gives the acoustics plane wave equation [Refs. 6,7] 
2 2 
2 p(x,t) aod p(x,t) 
a Spe 


2 
ot Ox 


(2.5) 


2. Spherical Wave Equaticn 
The next step in the process is to work backward from Eq. (2.5) and 
find an expression for the traveling wave in spherical coordinates. First, assume 
a point source with the pressure wave expanding radially. Instead of the single 
dimension x, Eq. (2.5) is three dimensional. Changing coordinate systems from 
linear to Cartesian and taking the Laplacian in place of the partial derivatives 
with respect to position Eq. (2.5) takes the form [Refs. 6,7] 


8 p(x,y,z.t) a 
ad V xy + Le 
: oY pQxyrzst), (2.6) 


ot 


. Assuming no losses, a uniform traveling wave is simply a function of 


distance r from the source and time, and is independent of angular displacement 














from the source. Therefore, it is logical to convert to a spherical coordinate 


system and use the spherical form of the Laplacian operator in Eq. (2.6). After 
some simplification [Ref. 6], Eq. (2.6) becomes 
2 2 
2 (ps) _ 29 GPG) 


; a (2.7) 
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where r is the radial distance from the source. This is the spherical form of the 
wave equation and is a function of only the distance from the origin, and time. It 
is the same as Eq. (2.5) with p(x,t) replaced by rp(r,t). 
Given a complex sinusoidal pressure function, a solution to Eq. (2.7) is 
given by 
p(r,t) = tm io) (2.8) 


where pmis the peak pressure amplitude at unit distance, and c is the nominal 


sound speed [Ref. 6]. 


B. RAYPATH DETERMINATION 

There are two fundamental reasons that an acoustic signal from a single 
source may have muitiple raypaths and corresponding arrival times at a single 
receiver. First, the ocean is a nonuniform medium. It varies in depth 
(pressure), salinity, and temperature. There are numerous currents, eddies and 
tidal effects. Second, it is limited vertically by a sea-air interface and a sea- 
bottom interface. Because of the first condition, sound speed is also nonuniform 
and is a functior. of salinity, temperature, and pressure. Fortunately, a great deal 
is known about the dependencies of sound speed in water and the behavior of 


sound waves at the sea-air and sea-bottom interfaces. Applying this knowledge, 


and Snell's Law, raypaths can be accurately determined. 








UT 


2. Sound. Speed 
Sound speed. formulation as a function of temperature, salinity, and 


pressure has been determined numerically, and is approximately [Ref. 7] 


c= 1449.2 + 4.6T - 0.055T“+ 0.00029T° + (1.34 - 0.01T)(S-35) + 1.58x10°P , (2.9) 


where c is sound speed, T is temperature, P is gauge pressure of a column of 
water, and S is salinity. For most applications salinity variations are small and 
can be neglected. Pressure is a linear function of depth. Temperature gradients 
and layers are found throughout the world’s oceans and are in general a function 
of depth. Warmer water tends to reside near the surface with colder 
temperatures at deeper depths. A certain amount of mixing can also occur to 
form isothermal and isosaline water. Temperature layers can also be formed due 
to currents and eddies. Sound waves can be thought of as tending to bend toward 
points of lower sound speed, as will become apparent in the next section on 
Snell's Law. 
2. Snell's Law 


A raypath through the ocean medium can be determined by considering 





the sound speed versus depth to be a continuous function, implying a 
continuously stratified medium. Changes in propagation path can be computed 
by treating the entire stratified medium as a set of n layers, and determining the 
refraction from medium | to medium 2, then 2 to 3, and so on to medium n. 
This is done through the use of Snell's Law, which is derived in several texts 


[Refs. 6,7]. Stated here, a sound wave (ray) will be refracted from medium to 


medium according to the relation 








sin(0,) sin(@,) sin(@,) 
Cy — C> — eae a : 


(2.10) 


where 0; is the incident angle, measured from the ray to the vertical, and Q3, ..., g 


6, are refraction angles. The constants c), C2, ... Cc, are the sound speeds in 





medium 1, medium 2, and medium n, respectively. It is assumed that the energy 
in the ray is constant through the boundary. 

As an example, consider the case where cj > cz. In this case, the sound. 
speed in-medium 1 is faster than that of medium 2. Rearranging Eq. (2.10) in 
the following manner 

Osin@y = sin(®.). (2.11) 
Cy , 
it can be seen that ratio of cz to c, is less than one. For the equality to hold, 02 
must be a smaller-angle than that of the incident angle, 6;. This means the ray is 
“bending” downward towards the medium with *he lower sound speed. Similarly 
the ray bends upwards for the case when c2>c.. In-other words, sound waves 
bend toward the region of slower sound speed. 

In reality, when sound impinges on a boundary between two mediums 
the sound waves are both reflected and transmitted. Application of Snell's Law 
remains the same for the transmitted wave, and the reflected wave has a 
reflection angle that is equal to the incident angle. The transmitted angle remains 
a function of the sound speeds of the mediuins on eithe: side of the boundary and 
is given by Eq. (2.10). However, the amplitude of ine transmitted wave dees not 
equal that of the incident wave, because con:...vation of energy must apply. 
Therefore, the amplitude of the reflected and iransmitted waves must sum to that 


of the incid-nt wave at the boundary. Assuming the density of medium 1 and 2 








to be the same, the ratio of the amplitude of the reflected wave to that of the 


incident wave is given by [Ref. 6] 
R C2sin(0,)— cc, sin(@,) 
ae a ae ARS (2.12) 
c2 sin(8,)-- .; #28.) 
where R is the reflected amplitude, 1 is th. iu° ¢ »"* amplitude, @; is the incident 
angle. @2 is-the transmitted angle, and cy, ai.d -. are-the sound speeds in medium 
1 and 2, respectively. Similarly, the za: + -€ the transmitted to incident 


amplitude is given by 


z Co Sin 8} 


ee 2.13 
Cy sin’ ,) +c; sin(@,)’ Gd> 


where T is the transmitted amplitude. For the case when the mediums do not 
have the same density, the sound speeds cj in Fq. (2.12) and Eq. (2.13) are 
simply replaced by 

Zj=Pj;C:» (2.14) 
where p,, and Z, are density and impedance of the-ith medium, respectively. 

For rays incident at ti. sea surface the sound wave is considered to be 
totally reflected, and the amount of reflection at the ocean bottom is based on the 
local bottom consistency. There are two special conditions that can be applied 
with Snell's Law. In the first case, when c; > c2 and the ratio of c; to cz becomes 
very large, the transmission angle @2 must approach zero. In this case, the wave 
transmits perpendicular to the boundary. In the second case, when cz > cj, there 
is a.condition when the transmitted angle must be negative to satisfy Eq. (2.10). 
However, 82 cannot be negative. This is the condition for total reflection and the 


reflected angle equals the incident angle, and 62 is zero. 





weet 





3. Travel Time 
Once the sound speed as-a function of depth has been determined, the 
ray path can be determined. Travel time can then be evaluated as an integral 


over the raypath. Travel time, Tj, for a particular ray i from point 2 to b is 


[Ref. 9] 


b 
_ ds 
T; -| C(x,y,2)’ (2.15) 


where sound speed is assumed to be independent of time t. If th assumption is 


‘made that 


C(x,y,Z) = Co(x,y,Z) + 5c(x,y,z) , (2.16) 


where Co(x,y,z)-3s an assumed background sound speed field, and 5c(x,y,z) is an 


unknown sound speed field, Eq. (2.15) can then be written as 


b 
b 
d dc(x,y,z) ds 
Tj= To + 8T;= [ CRT 2 ae (2.17) 


2 
7 J, Co(X,Y,2) 


for dc << co. The travel time T; has been broken down into a travel time To, cue 
solely to the assumed sound speed field co, and 6T;, a linear function of the 
unknown sound speed field perturbation dc(x,y,z). Variation in the arrival time 
T, is a j-erturbation of the arrival time due to perturbations in the sound speed 
along the path By measuring these variations and using linear inverse 
mathematical techniques it is possivie tc determine some desired characteristic of 


this geophysical problem [Refs. 2,9]. 

















III. M-SEQUENCES AND HADAMARD PROCESSING 


A. INTRODUCTION 

An important parameter in undervater acoustic tomography is the arrival 
time of a transmitted signal. A very usefu: means for measuring arrival time is 
to apply an impulsive excitation to the system, in this case the ocean 
‘environment. An impuise can be generated by explosive or implosive sources 
cheaply and easily, but suffer from unevenness in the frequency spectrum, and 
repeatability. Another method ‘avolves the use of pseudorandom noise. 
Psuedorandom noise is deterministi atu is formed by a sequence of ones and 
zeros known as a maximal-length sequence or m-sequence [Refs. 10,11]. M- 
sequences are typically transmitted by using the ones and zeros of the sequence 
({1,0] mapped to [-1,1]) to phase encode a carrier signal ['%efs. 12,13,14]. 

M-sequences are generated using a simple binary shift register that is formed 
«cording to the desired primitive polynomial. Since the signal being transmitted 
is known, a simple autocorrelation can be formed uring a bank of matched 
filters. The resulting correlation is impulse-like, and has a shorter duration than 
the originally transmitted signal and has a much larger amplitude, which makes 
measurement of arrival time an easier task. 

The following sections give a brief description of shift register fundamentals, 


m-sequences, and fast Hadamard processing of m-sequences. 














B. SHIFT REGISTER FUNDAMENTALS 

A general sequence shift register generator is shown in Figure 3.1. A 
sequence of ones and zeros can be output from any one of the registers starting at 
any register state, except zero. In what follows all output sequences are taken 
from the Oth register. The sequence of bits generated is periodic and will repeat 
with some period. L, based on the structure of the generator, as described by a 
polynomial such as 

g(D) = b,D"+b,, ,D"'+...+b;D+by» (3.1) 

where D is the unit delay and by are the feedback weighting coefficients. The 


weighting coefficients take on values of one or zero indicating connection or no 
connection to the kth register. All arithmetic operations for polynomials are 
performed using finite-field arithmetic (modulo two) and are described in more 


detail in Reference 10. 





Figure 3.1: General sequence shift register generator for Eq. (3.1). 
An example of a third order sequence shift register generator is shown 
in Figure 3.2, and its corresponding generating polynomial is 
e(D)=D3+D+1. (3.2) 
In this case Eq. (3.2) happens to be a primitive polynomial, where a primitive 
polynomial is any polynomial that will not repeat its register state until after 


2.1 delays, where n is the number of delays in the shift register generator. 














Therefore, a shift register that is defined by a primitive polynomial will generate 
a sequence that is of maximal-length and is known as a maximal-length sequence 


or m-sequence. 





Figure 3.2: Shift register realization of Eq. (3.2). 
Continuing with this example, the sequential output from the registers of 
s igure 3.2 is then given as in Figure 3.3, where the initial register contents is 
arbitrarily set to ay=1, a;=0, ag=0 [Ref. 9]. 
The corresponding m-sequence is then obtained by taking any one column 
from the register states top to bottom (or bottom to top). Note, that the 
combination of [000] never occurs. This is because an initial value of zero will 


not allow any transition in state and can-be easily verified by inspection. 











Cycle a, a a 


1 0 
1 1 0 0 
2 0 1 0 
3 0 0 1 
4 0 
5 1 1 1 
6 1 1 0 
7 0 1 ] 
8 1 0 0 


Figure 3.3: Shift register contents when generating a third order m- 
sequence. The eighth cycle shows that the register begins to repeat. 


The characteristics of the m-sequence is unchanged whether it is transmitted 
in the forward or reverse direction, but when performing the fast Hadamard 
Transform, reviewed in the next section, the direction in which the sequence is 
transmitted is significant. Therefore, the top to bottom sequence will be 
designated the "forward" code and the bottom to top the "reverse" code, 

forward m= 0011101 
reverse m= 1011100. (3.3) 

M-sequences have several desirable characteristics. Its correlation function 
is triangular in shape, see Figure 3.4, and short in duration. It is deterministic, 
once the polynomial is determined all output is known. Jt can be implemented 
easily and is periodic. However, one major drawback is that its autocorrelation 


requires N multiplies, and since the arrival time is never known, the input signal 


must be correlated with all N shifted versions of the original sequence which 











rt 


requires N* multiplies. This can be overcome by use of the fast Hadamard 


transform (FHT), discussed in the next section. 








Figure 3.4: Autocorrelation function for a m-sequence of Jength N. 


C. THE FAST HADAMARD TRANSFORM 
The autocorrelation of the input data sequence with all possible shifted 


versions of the original m-sequence can be written as 


. O01i10ilfa 
40014410] ]>d 
mp- (0200424]/¢ 
=14010011 2 
4140400il/e G4) 
44i0100]/f 
0441010] ]/g 


where M is a matrix whose rows are the shifted versions of the forward code and 
D is the input data vector. This product requires N* multiplications. To reduce 


this number, a Hadamard transform can be used. The third order Hadamard 


matrix is given by 








oi: Gite a ee a a Ce 
4-1 4-1 1-4 4-4 
1 1-1-1 1 4-4-1 

Ha | trt-1 4 1-4-1 4 
41441 1-4-1-1-1 (3.5) 
4-1 4-1-1 1-1 1 
-4 1-4-41-4-1 44 
4-1-1 1-1 141-1] > 

where H can. be formed recursively from 
Hy= [1] ; Hie = | Hi Hi (3.6) 
i7Hy] 


and by performing a simple mapping of (1,-1) to (1,0). H can be written in a 


form that is easier for most people to work with [Refs. 15,16,17], 
00 0 0 


(3.7), 


cCOoCDOOCOD 
POPORPOR 
PPOOPROO 
OPPOORP 
PRERPCOOO 
DRPORRPORO 
CORRRPPRPOD 
FPOOPRPOP 


It is can. also be shown {Refs. 16,17] that the matrix H can be factored into 


the product of two matrices consisting of a binary count, in this case from 0 to 7, 


| G8) 


PPRPRPOOOO 

PRPOORPRPOO 

FPOPORORPS 

Levens 

oro 
ERO 
oOF 
POR 
ORE 
PPE 


The matrix M, Eq. (3.4), can also be factored into two matrices, L and S, given 


by [Refs. 16,17] 
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It is easily verified that 





M=LS. (3.10) 


Note that by adding a leading column of zeros to S, forming S', and a leading 
row of zeros to L, forming L'. All possible combinations of ones and zeros are 
formed in L’, S’, as in A, with the differences between the matrices being the 
order in which they occur. It is straight forward to find a matrix P such that 
S'=A7 P and another matrix U such that L'=UA. Combining these results maps 
M' to the Hadamard matrix as [Ref. 16,17] 

M' = L'S' = UAA™P = UHP, (3.11) 
where M' is the M matrix with an appended column and row of zeros in the first 
row and column. Recall that the correlation was performed by forming the 
product of the M matrix and the data vector D. By replacing M with M' and 


forming a new data vector 


a 
b 
c 
D' = d 
5 (3.12) 
f ) 
g 
| 
and after substituting Eq. (3.11) for M’ the correlation becomes 
R’ = M'D' = UHPD". (3.13) 
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From this it can be seen that the matrix M is not needed and the Hadamard 
matrix can-be used in its place. By using the Hadamard matrix in the form of 
Eq. (3.5) and forming a product with D', the result becomes a sum of the 
individual terms of D' with + and - weights and has the form 
| atbt+tctdtet+ft+gth 
a-b+c-dt+e-f+g-h 
atb-c-d+e+f-g-h 
ad Pane a aA (3.14) 
a-b+c-d-et+f-gth ; 
atb-c-d-e-ftgth 
a-b-ctd-e+ft+g-h 
which when written in the form of a flow graph having the same form as the 
well known fast Fourier transform shown in Figure 3.5, where the complex 
multiplies are set to one (Refs. 9,16,17]. This is known as the fast Hadamard 
transform (FHT). 

All that is left is to determine the permutation matrices P and U. By looking 
at how P and U must be formed, it can be seen that P must have ones at the 
indices [Ref. 15,16] 

ROW 0 1 2 3 4 5 6 7 

COLUMN 0 2 1 4 6 t 3 a3 
and U must have ones at the indices 

ROW 0 4 2 1 6 3 7 5 

COLUMN 0 1 Z 3 4 5 6 Ts 
which correspond to the binary values of S' and L’. It turns out that the matrices 
S' and L' provide all of the necessary information. All that is needed to perform 


the desired multiplication is to permute the input data vector according to S', 


form the FHT, and permute again according to L’. L and S are generated from 
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the primitive polynomial that defines the m-sequence using the generators shown 
in Figure 3.6 for the forward code and Figure 3.7 for the reverse code [Ref. 15]. 
The modification to form S' and L' simply involves adding leading zeros to S 
and L. Since the correlation is always over 2"-1 data points, a leading zero is 
added to the input data vector D forming D' adding no new data. After 
processing, D' holds the correlation function with the zero position containing 
the average level, which may be used to remove the bias in the autocoueiaiion 


function [Ref. 11]. 


Basic Element 


A A+B 


B A-B 





a+b+c+d+e+f+g+h 
a-b+c-d+e-f+g-h 
a+b-c-d+e+f-g-h 
a-b-c+d+e-f-g+h 
a+b+c+d-e-f-g-h 
a-b+c-d-e+f-g+h 


a+b-c-d-e-f+g+h 





a-b-c+d-e+f+g-h 


Figure 3.5: Basic fast Hadamard transform element for cascading 
additions and the full diagram for an eight point FHT. 


21 




















Column i Column 2 Column n 


Row 1 





(b) 


Figure 3.6: (a) L generator and (b) S generator for the forward 
code. 


22 











Row 4 


Row 2 


Row n 





(b) 


Figure 3.7: (a) L generator and (b) S generator for the reverse code. 














IV. TIME-VARYING DOPPLER PROCESSING 


A. PHASE ENCODING OF M-SEQUENCES 
A simple method in which m-sequences can be transmitted is to phase encode 
the digits using a mapping of (0,1) to (1,-1) and appropriately selecting a phase 
angle which will maximize the signal-to-noise performance. The phase angle 
that optimizes this is 
y= tan '({N), (4.1) 
where y is the phase modulation angle, and N is the length of the sequence [Ref. 
11). However, for the Heard Island Experiment, y = 45 degrees is used to 
simplify SNR estimation. The transmitted signal s(t) takes the form 
s(t) = Acos(2nf,t + M(t)y), (4.2) 
where A is the amplitude, f, is the carrier frequency, and M(t) takes on values of 
+1 or -1 depending on the m-sequence code. The minimum length of time M(t) 
remains constant is the digit duration d. 
It is normal practice to chose an integer number of cycles, Q, of the carrier 


frequency per digit. 


B. SIGNAL PROCESSING WITH ZERO DOPPLER 


At the receiver the signal is normally sampled at an integer multiple of the 


carrier frequency f,, i.e. with fy = mf,, where m is chosen to be an integer 
greater than two, to ensure that the Nyquist criterion is satisfied, and f, is the 
sampling frequency. Sampling in this fashion simplifies processing and 
interleaving of data, since there is an integer number of data points per digit d. 


The received signal is r(t) = s(t - T,), assuming no attenuation, no dispersion, and 
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a single raypath. The sampled input r(n) is then den odulated to remove the 
carrier from the signal. The received signal r(n) will arrive at some unknown 
arrival time Tj, and must be multiplied by both the sine and cosine functions to 
form the in-phase and quadrature components p(n) and qin), respectively [Ref. 


18]. For p(n) this then gives 











p(n) = Aco + me - T)yoeos— 
p(n) = . [cos(M- -T)y) + cos * MG - Typ), (4.3) 
and similarly for q(n) 
a(n) = 5 [sing - Ty) - sine +MG-- TDW) (4.4) 


The high frequency components are then removed by lowpass filtering with the 


cutoff frequency greater than the digit bandwidth. The resulting waveforms are 


pln) =F cos(M(q-- TW) 


q(n) = as sinMG-- T)y), (4.5) 


which are constant over the length of the digit of M(n). 

Because of the lower digit bandwidth, decimation in time may be performed 
without any loss of information [Ref. 19]. This significantly reduces the data 
rate and storage requirements of the processed data. Further processing follows 
the scheme of the previous chapter. The input data vector is permuted, the FHT 
is performed, and permuted again to the correct order. Multiple samples per 
digit must be accounted for by interleaving when performing the FHT on a given 


sequence length. Finally, magnitude and phase are computed in the normal 


fashion. 




















The autocorrelation function, see Figure 3.4, for m-sequences has a negative 
DC level for all digit positions.except the first. A correction may be made which 
will adjust this level-to zero. Knowing the phase modulation angle, the DC bias 
-can be removed by adding the correction 
Coor = filave ? (4.6) 
where Cor is the complex level adjustment, Lay, is the average DC level from 
the FHT, and f, is given by the relation [Ref. 11] 
6 = ED tanty) +N - tan’(y)) | 
N* + tan“(y) (4.7) 
where N is the length of the m-sequence and y is the phase modulation angle. 
Since all parameters of Eq. (4.7) are known, f, is computed once at system 
initialization and Coo, requires one complex multiplication for each data length 
processed. 
A block diagram showing the process is-given in Figure 4.1. Note that in 
this case a Butterworth filter of order ten was used in the passband and a 
Chebychev filter of order five was used for the lowpass filter. Any filter type 


can logically be substituted. 
















10** order | 
Buttervor thi 
BP Filter. 
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cos eC @- sin(2nf,n/f.) 
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Magnitude 


Figure 4.1: Block diagram of signal processing with zero Doppler. 











C. SIGNAL. PROCESSING WITH DOPPLER 

Generally in acoustic tomography the receiving or transmitting ship 
must maintain way and/or buoys are acted on by the forces of currents, tides, 
wind etc., all of which cause relative motion between the receiver and 
transmitter. The Doppler shift that is then imparted on the signal must be 
corrected for at the receiver. For the relatively slow and uniform speeds 
encountered in acoustic tomography it is sufficient to compensate for first order 
Doppler, and neglect acceleration effects [Refs. 20,21]. 


When Doppler is present a transmitted_signal s(t) may be received as 
Vv 
1(t) = s((1+=)t- Til, (4.8) 


where v is the velocity in meters/second, c is the nominal sound speed in water in 
meters/second, and Tj is the delay from transmitter to receiver as in section B, 
where attenuation, dispersion, and a single raypath are assumed. The magnitude 


of the spectrum of r(t) can be given as 


ROIs E 





] 
GY) t (4.9) 





Sampling at an integer multiple of the carrier frequency assuming zero Doppler 


and replacing the time dependency with the discrete index n gives 


r(n) = -T)), (4,10) 


With Doppler, the received signal is then given by 














n aE 
r(n) = s[(— ) - Ty), (4,11) 
Adjusting the sampling frequency by the same Doppler shift gives a new 


sampling frequency of 
f=f,(1+V/), 


(4.12) 
which forms a new received signal 
n . + 2, 
r'(n) = s{[(—=—)- T] 
fs (4.13) 
Substituting Eq. (4.12) into Eq. (4.13) yields 
r(n) = s( -T;) =r), 
fs (4.14) 


At this point the results are identical and processing of the data can proceed as in 
the zero Doppler case [Ref. 20]. 

Two methods can be employed to perform the adjustment to the 
sampling frequency as discussed. The first is to use a bank of sampling devices 


that sample at a set of sampling frequencies that cover the desired Doppler range 


such that each device samples at frequencies f,,,f,o,...£,. This method requires a 


great deal of hardware and is inflexible and expensive to implement. The second 
method is to frequency shift the received signal and interpolate between samples 
at the new sampling frequency to predict the value of the signal as if the desired 
sampling frequency had been used [Refs. 20,21]. 

Doppler resolution is 1/T, Hz, where T, is the analysis interval] [Ref. 
20]. For example, a N=255 length m-sequence with Q=5 and f,=57 Hz giving a 
T, of 22.368 seconds and a resolution of 0.04471 Hz, which from Eq. (4.12) 


corresponds to a Doppler shift of + or - 2.3 knots. Figure 4.2 is an ambiguity 


plot derived from this example. The plot shows signal arrival time for the given 

















m-sequence versus signal Doppler over-a + or - 5 knot range. As can be seen in 
the plot signal strength decreases as the carrier frequency moves away from the 
zero Doppler case and is not detectable at about 2 knots on either side of center. 
Frequency shifting in the frequency domain is equivalent to multiplying 
the samples or the demodulates in the time domain by exp(j2mnf,/f,) [Ref. 20], 
where f, isthe Doppler shift of the carrier signal in Hertz. The Doppler shift f, 


is-determined directly from the expected Doppler and is given by 


Vv 
f= che (4.15) 
For a sampling frequency that is four times the carrier the multiplication 


exponential can then be related to the Doppler speed, v (meters/second), and the 


multiplication of the demodulates in the time domain is by exp(jznv/2c). The 


ambisuity Ploc 
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22.368 sec 
Figure 4.2: Ambiguity plot of Arrival Time vs. Doppler for a 


theoretical sequence with N=255, Q=5, and fe=57 Hz. 














data is then interpolated at the new sampling frequency corresponding to the 
original sampling frequency shifted by the same.amount of Doppler as given by 
Eq. (4.12). 

The Doppler processing scheme in this thesis is performed as follows. 
Frequency shifting is performed by combining the shift with the demodulation 
process. This reduces the number of complex multiplies by combining the 
modulation angle and the frequency shift via simple addition, see Figure 4.3. 
The resampled signal is then linearly interpolated [Ref. 21]. Since the Doppler 
shift is unknown, except within some range, each data set is processed over the 
entire Doppler range, where the Doppler step size is set in knots. Frequency 
shifting and then interpolating results in the following expression for the 


resampled signal [Ref. 21] 


a jarfat, Me. j2nfaty jarfat, 
s(t;+h) = s(t,)e + ee [s(tye s(t,Je iF (4.16) 


where t,< t,+h<t, and time replaces the discrete index n for clarity. The times t, 
and t, correspond to sampled data points, and h is the time spacing from t, to the 
point to be interpolated. The position h, varies from interpolation point to 
interpolation point, according to the sampling period I/f,', as time progresses. 
Other techniques for interpolating between samples are available but are slower 
and more complex and will not be discussed here. 

At this point processing continues as in the zero Doppler case. Figure 
4.3 shows the process in block diagram form with an additional block for 
performing coherent averaging to increase signal processing gain, if required or 
desired (Ref. 18]. 


Figure 4.4 shows the Doppler resolution for a transmitted signal with 


-1.4 knots of Doppler, using the parameters given in the above example without 











noise, and. was processed by this scheme. Note that these results are plotted on a 
different scaling than that of Figure 4.2 but show the same results. The Doppler 


resolution is again approximately + or - 2 knots as expected. 
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Figure 4.3: Block diagram of signal processing with Doppler and 


coherent averaging incorporated. 
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Figure 4.4: Ambiguity plot for Arrival Time vs. Doppler for 
sequence length N=255, Q=5, and fe=57 Hz, processed by SEQREM 


program with no signal noise. 








V. RESULTS AND CONCLUSIONS 


A. RESULTS 

The objective of this thesis was to develop a program in the C programming 
language that can process any m-sequence at any carrier frequency over any 
acceptable Doppler range. This objective has been fully met. The package 
SEQREM has been developed. which will scan over any specified Doppler range 
and specified Doppler step size. In cases where no Doppler is expected the 
program skips the Doppler scanning procedure, reducing run time. 
Programming allows for any carrier frequency, any type filter of any order, up 
to twenty, and any size m-sequence, with dynamic memory space allocation for 
array manipulation, limited only by the available computer memory. It permits 
coherent averaging of sequence segments as specified at initialization, improving 
processing gain, and computes and removes the DC bias in the autocorrelation 
function. Decimation in time is performed which significantly reduces off-line 
storage requirements for processed data. 

Appendix A contains a complete set of processing results for all three cases, 
closing, opening, and zero Doppler. In all of the plots a 255 digit sequence is 
used with Q = 5, and a carrier of 57 Hz, corresponding to one of tie signals to 
be transmitted in the Heard Island Experiment. The sampling frequency ‘s four 
times the carrier frequency or 228 Hz. The simulated data assumes white 
Gaussian noise. The plots are labeled with the Doppler that was processed and 


with the sequence period, adjusted for the expansion or compression, according 














to the Doppler bin. Decimation in time was performed at a ratio of ten to one. 
Peaks correspond to the center of detection of the signal arrival within the 
period. The leading edge of the peak is the actual arrival time and can be 
extracted by an edge detection program. Sequential frames are aligned one 
behind the other and represent increasing time along the y direction. All plots 
are in magnitude. Phase plots are omitted for brevity. Some specific results are 
contained here. 

Figure 5.1 shows the detection of a zero Doppler signal that has been shifted 
in time and has some arbitrary arrival phase. The signal is shifted approximately 
three digits with an SNR = 0 dB. 

The sequence removal process is linear and can detect multiple signal 
arrivals. Figure 5.2 shows this clearly for two arrivals spaced 27 digits apart 
(approximately 6 seconds) with SNR = -15 dB. Figure 5.3 shows the same signal 
with one digit separation, and the signal arrival times are still resolvable. 

The next figure, Figure 5.4, demonstrates the processing gain from coherent 
averaging. The data is the same as that in Figure 5.1, but with each segment 
output averaged over three frames, and a larger shift in arrival time. Note the 
increased SNR as compared with Figure 5.1. 

Figures 5.5 and 5.6 show the detection of two Doppler shifted signals. The 
first is for an opening situation with -1.4 knots of Doppler and -12 dB SNR. The 
second is a closing situation with 3.4 knots Doppler and a SNR of -15 dB. 


Graphs of all Doppler bins searched are include in Appendix A. Appendix B 


gives a brief description of the program SEQREM. 




















ARRIVAL TIME PLOT 
@.0@ KNOTS DOPPLER 





Figure 5.1: Zero Doppler signal. N=255, Q=5, fe=57 Hz, 0 dB SNR. 
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Figure 5.2: Zero Doppler signals, 27 digits apart. N=255, Q=5, 
fe=57 Hz, -15 dB SNR. 
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Figure 5.3: Zero Doppler signals, 1 digit apart. N=255, Q=5, fe=57 
Hz, 0 dB SNR. 
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Figure 5.4: Zero Doppler signal. N=255, Q=5, fe=57 Hz, -15 dB 


SNR coherently averaged over 3 frames. 
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Figure 5.5: -1.4 knot Doppler signal. N=255, Q=5, fe=57 Hz, -12 dB 
SNR. 
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Figure 5.6: 3.4 knot Doppler signal. N=255, Q=5, fe=57 Hz, -15 dB 
SNR. 
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B. CONCLUSIONS 

The results of the previous section show that the software developed for this 
thesis operates as desired. Sequence detection is performed as predicted, via 
software compensation of sampling frequency, for any reasonable Doppler using 
linear interpolation techniques. Linear interpolation may not be the optimum 
interpolation method, but it has been demonstrated that it does work and is 
simple and faster than more complex techniques. This software will enable the 
Naval Postgraduate School (NPS) to participate in the Heard Island Experiment 
and will provide an important processing capability for any future experiment, 
and most importantly, a Doppler capability that was not previously possessed by 
the Underwater Acoustics Systems Laboratory of the Electrical and Computer 


Engineering Department at NPS. 


C. ADDITIONAL WORK 

Despite all of the advantages of this software, further work is necessary. 
The program as it stands, operates close to real time on individual segments, 
given the inherently low data rates used in acoustic tomography, but resides ona 
ULTRIX-32 V2.0 operating system with a VAX 11/785 processor. This is a 
time sharing system and is not designed for real time processing. Because of 
this, long data sets end up with a reduced priority as time progresses, causing 
slower execution from the user's perspective. 

This problem will be remedied with the delivery of a Concurrent Computer 
Corporation VME Native-Mode Data Acquisition System in November of 1990. 


The intent is to transfer this software from the ULTRIX-32 system to this 


machine. The VME supports multiple processors and can process 14 MIPs at 33 











MHz. It uses a memory cached system of up to 128 MB. It can also support up 
to 7 separate graphics heads. This system is capable of real time operation with a 
UNIX type operating system and will be ideal for use with this programming 


package. Given this-system, the work that needs to be done includes: 


1. Transfer of the existing program and related functions to the VME 
system. This includes the SEQREM program, initialization 
functions, FHT function, and magnitude computation routines. 


2: Modification of the programs to incorporate real time processing in 
a parallel vice serial-structure. Data input to the system is first 
read, processed, stored to file, and then the next segment is read 
and so on. A parallel scheme is needed to allow reading of the next 
data set while computation is progressing on the current while the 
results of the last set are being displayed. 


a: Design of an I/O system to handle data input and output channel(s). 
Without the VME available, it was necessary to assume input data is 
maintained in text file. Output is also to a file with sequential 
Doppler searches appended to the same data file. The software also 
assumes single channel input. It would be desirable to handle 
multiple channels. 


4. Development of a display software/system to present the data in an 
convenient form, such as a waterfall type display. Currently, 
plotting of data is performed off-line, after all computations are 
completed. There is little flexibility for plotting various size data 
sets and real time display is not possible. 


3 Altering the Doppler bin search from a sequential to a parallel 
operation. Currently in the data processing section, Doppler 
processing is performed in increments via a standard "for" loop. 











APPENDIX A 


The results of Doppler processing for the three possible cases are given in 
the following figures. Figure A.1 shows a zero Doppler signal. No scan is 
made, since no motion is expected. Figure A.2 is a closing Doppler situation. 
Doppler search is over a + or - 3.5 knot range, in 0.5 knot steps. Figure A.3 is 
the opening case with a search over + or - 3.0 knots, also in 0.5 knot increments. 
In each figure the following parameters apply: primitive polynomial = 537s, N = 
255, carrier = 57 Hz, w = 45 degrees, initial state = 13, Q=5. White Gaussian 
noise is assumed. Plots are displayed as arrival time versus time. This data 


simulates a signal to be transmitted in the Heard Island Experiment. 
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Figure A.1: Zero Doppler signal. SNR=0 dB. 
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(b) 
Figure A.2 : 3.4 knot Doppler signal, 0.5 knot steps. SNR= -15 dB. 
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Figure A.2 : (cont.) 
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: (cont.) 
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Figure A.2 : (cont.) 
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Figure A.3 : (cont.) 
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APPENDIX B 


The following programs make up the SEQREM program or are associated utility 
programs. Since sufficient documentation is included with each program, only a very 
brief description is provided for each. This software may be obtained on floppy disk 
by contacting Professor James H. Miller at (408) 646-2384 or by writing to his address 
in the initial distribution list on-page 98. 


A. MAIN PROGRAM 
1. MACROFILE 


This is a a short file that holds definition macros and global variables. 


/* MACROFILE includes the necessary libraries and declares global constants 
and global variables for general use. 


Variables: 
scram ~ pointer to indices to scramble input data vector prior to FHT. 
unscram ~ pointer to restore data vector after FHT performed. 
law - polynomial law used to generate M-sequence. 
degree - the degree of law. 
initiai - initial register load for generation of the sequence. */ 
finclude <stdio.h> 
#include <math.h> 
#include <malloc.h> 
#define BTRMAX 11 
#define PI 3.1415926536 
#define FALSE (unsigned) 0 
fdefine TRUE (unsigned) ! FALSE 


unsigned *scram, *unscram, law, degree, initial; 
char *malloc(); 


2. SEQREM 


This program is the root of all the processing functions for sequence removal. 


It uses standard file J/O and may be used with any standard input numerical format. 
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The program call is SEQREM, and the user is prompted for the following 
information before processing may begin (Note: some parameters are requested from 
within initialization functions, however, this is transparent to the user): 

1. Data file containing formatted filter coefficients. 
2. Primitive polynomial defining the m-sequence. 
3. Initial register load. 
4, Transmission direction, i.e. forward or reverse code. 
5. Carrier frequency. 
6. Sampling frequency (checks for four times carrier). 
7. Phase modulation angle. 
8. Cycles per digit, Q. 
9. Desired time decimation, if any. 
10 Coherent averaging desired, if any. 
11. Expected Doppler range (knots). 
12. Doppler search step size (knots). 
13. Input data file. 
14. Output data file. 
Assumptions: 
1. Sampling at four times carrier. 
2. One period of m-sequence processed at a time. 
3. M-sequence is taken from the least significant register. 


4. Decimation ratio is a rational. 


/* SEQ_REM performs sequence removal of a phase encoded M-Sequence. t is designed to work 
for any maximal length sequence as defined by a primitive polynomial. The polynomial, 
carrier frequency, sampling frequency, etc. are input by the user. All filtering and 
demodulation is performed, but filter coefficients must be stored in an externa. file in the 
format described by the function GET FLT_COEF(). For ease of use, these may be generated off 
line by the MATLAB program SET_FILTER.M. 














M-Sequences my be transmitted in the forward or reverse direction. 


? - is used as a wild card. for brevity and may indicate 1,2,a,b etc. 
in the following variable and file descriptions. 


Variables: 
c,j},1,k,n,count - integer counters used in various loops. 
dir,yesno - variables used for user inqueries. 


num_pts - number of points processed per digit, before and after 
decimation. 


num_coh_avg ~ number of frames to be coherently averaged. 
vect_len - numoer of points in a sequence before decimation in time. 


vect_len_inter - number of points need to be read for interpolation 
of one segment length at the new sampling freq. 


step - step size determined for decimation in time. 
cycles - number of carrier cycles per digit, must be integer number. 


first_run - flag to determine if a segment is the first one to be 
processed in the current doppler bin. 


end ~- flag to determine when end of file has been reached. 
seq_len - number of digits in a M~sequence. 
highpass_?, lowpass?_? - pointer to storage for filter coefficients. 
indatal,indata2 ~ array for data storage and processing storage before 
decimation in time. Indatal is initially input 
and then is imaginary part and indata2 is real part 


{after demodulation). 


h_data_im,h_data_re ~- arrays for real and imaginary vectors for storing 
scrambled data for FHT processing. 


re data,im_data - arrays for real and imaginary data storage after data 
unscrambling. 


avg_re, avg_im - arrays for real and imaginary data storage when 
coherent averaging is performed. 


mag, phase - array for magnitude and phase of processed data. 
fc,fs - pointers for carrier frequency and sampling frequency. 


delvl_im,dclvl_re, pdstal_re,pdstal_ im - used to remove correlation 
de bias. 


ifact,rfact,ph_ang ~- factors and phase angle used to compute 
correction to de bias based on phase modulation 
angle. 


den - tmeporary storage used in above calculations of dc bias. 
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snd_vel - input average sound velocity to be used. 





{s,Ta,Ti - sampling period, sequence period, and interpolated sampling 
period in doppler interpolation. 


- time _s,time_i ~- current relative sample time and interpolation time, 
which are used to compute h. 


h - distance from tire a to sample to be interpolated where 
“ time a < h time b. 


ptri,ptr2 - used to save position in array space when performing 
doppler search. 


doppler - current doppler shift being processed (m/s). 


doppler _start,doppler_end ~ beginning and ending of doppler search 
interval (m/s). 


doppler_step,doppler bin - doppler step size used in scan, and doppler 
converted to xadians for demodulation. 


Functions: 


init_had() - gets relevant parameters for computing scambling and 
unscrambling indices. 


fwd_had(),rev_had() - compute indices for foward and reverse trans- 
mitted sequences. Both return polynomial 
degree from input. 


lowpass?(),highpass() - identical functions for rerforming filtering 
“ operations. 


get_fit_coef() - retrieves filter coefficient data from user provided 
file. 


demodulate() - performs demodulation for input signal to dc. Has an 
argument doppler bin for doppler frequency shifts. i 


} 


hadamard() - performs the FHT on scrambled data. Identical to F 
with multiplication factors omitted. 


1 
»! 





mag_phase() - computes the magnitude and phase of the processed cata.*/ 


Zinclude "macrofile.h" 


main () 
{ 
int c,j,k,n,i,dir, cycles, num_pts, step, vect_len, vect_len_inter,count, num_coh_avg; 
unsigned init_had(), rev_had(),secq_len, first_run, end; 
void get_flt_coef(),hipass(),lowpassi{) ,demodulate() ,hadamard(),mag_phase{); 
void lowpass2 (),lowpass3(); 
float *hipass_a,*hipass_b, *lowpass_a, *lowpass_b, *lowpassi_a, *lowpassi_b; 
float *indata2, *h_data_re, *h_data_im, *re_data, *im_data, *fs, *fc; 
float *mag, *phase, *avg_re, *avg_in,dclvl_re,dclvl_in,pdstal_re, pestal_in; 
float *data_i re, *data_i_im, rfact, ifact, den, ph_ang, doppler, sad_vel; 
float doppler_start,time_s, time_i, h, Ts, Ta, Ti, *temp_ptri, *temp_ptr2; 

- float doppler _end,doppler_bin, doppler_step, *indatal, *new(); 
int yesno; 
















































char infile(21), outfile[21); 
FILE *fpl,*fp2; 


/* Initialize filter coefficients for a low and high pass filter. */ 
/* These values are computed and stored in a file off line. x/ . 


count = 0; 


/* Allocate memory Vor fc, fs, and filter coefficients to be used. 10th 
order BUTTERWO4..! bandpass, and 5th order CHEBYCHEV lowpass filters 
must be used. Filter coefficients are retrieved and filters are 
initialized. Memory allocated to filter coefficients is freed after 
initialization. */ 


fc = new(1i); 

fs = new(1); 

-hipass_a = new(BTRMAX) ; 

hipass_b new (BTRMAX) ; 

lowpass _a = new(BIRMAX); 

lowpass _b = new(BIPMAX); 

lowpassi_a = new(BTRMAX) ; 

lowpassl_b = new(BIRMAX) ; 

get_flt_coef (hipass_a,hipass_b, lowpass_a,lowpass_b, lowpassi_a, lowpassi_b); 
hipass (hipass_b, hipass_a, BTRMAX) ; 
lowpassi(lowpass_b, lowpass_a, BTRMAX-5) ; 
lowpass2 {lowpass_b, lowpass_a, BTRMAX-5) ; 
lowpass3 (lowpassi_b, lowpassi_a, BIRMAX) ; 
free (hipass_a); 

free (hipass_b); 

free (lowpass_a); 

free (lowpass _b); 

free (lowpasy1_a); 

free (lowpassl_b); 


/* Set sound velocity to desired value, for given conditions. */ 
print’ ("\nEnter the Average Sound Velocity (meters/sec.)."); 


printf ("\n\n Sound Velocity: "); 
scant ("%i", &snd_vel) ; 





/* Injtialize parameters for M-Sequence, and set transmit direction. 
Set seurambling and urscrambling indices accordingly. */ 


initial = init_had(); 
prinef("\nIs the M-Sequence transmitted in the foward or reverse\n"); 
printf ("direction? \n"); 
nrintfé ("\n (forward=0/reverse=1j): "); 
scanf ("sd", &dir) ; 
1f£ (dir == C) 
degree ~ fad_had(law, initial, scram, unscram); 
else 
degree= rev_had(law, initial, scram, unscram) ; 
seq_len = (l<<degree)~-1; 


/* Compute offset to remove DC bias from correlation. */ 
printé("\ninter the Phase Modulation Angle used to encode the signal.\n "); 


print£ ("\n Phase Angle: ")s 
scanf("%f£", Sph_ang); 
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ph_ang = ph_ang*PI/180; 

den = seq_len*seq_len + tan(ph_ang) *tan(ph_ana); 
ifact = (seq_lent1) *tan(ph_ang) /den; 

rfact = -(seq_len-tan(ph_ang) *tan(ph_ang)) /den; 


. /* Input of transmission parameters and demodulator is initialized. */ 





printf ("\nEnter the Carrier Frequency used.\n "); 
printf£("\n Fe: "); 
- scant ("SE",&fc(0)); 
print£("\nEnter the Sampling Frequency. Sampling Frequency must be an"); 
printf ("\nfour times the Carrier Frequency.\n"); 
printf ("\n Fs: "); 
scanf ("%£", &£s[0)); 
1£ ( *fc/(*fs) != 0.25) 
( 
printf ("\nSampling Frequency will not work with this program. Exiting!\n"); 
exit (1); 
} 
printf ("\nEnter the number of carrier cycles per digit. An "); 
printf ("\ninteger number of cycles must be used.\n"); 
printf ("\n Q: "); 
scanf ("%d", &cycles) ; 


/* Compute sequence length in seconds (Ta) and sampling period. The number 
of points generated per digit is computed and provided as an aid in 
determining the decimation in time to be used, if desired. The total 
number of points per sequence period is computed (vect_len).*/ 


num_pts = (*fs/(*fc))*cycles; 
vect_len = num_pts*seq_ len; 
Ta = (float) (seq len * cycles / *fc);/* computes the sequence length in sec.*/ 
. Ts = 1/(*fs); /* compute sampling period. */ 
printf("\nThere are td pts per digit. For Processing Savings ",num_pts); 
printf ("\ndecimation in time is performed. The default used is one point"); 
priitf("\nfor each cycle in the digit. In this case %d points are",cycles) ; 
printf ("\nprocessed per digit. (i.e. every %d th point.)", (num_pts/cycles) ); 
printf ("\n\n Use default? (yes=1/no=0): "); 
scanf ("Sd", &yesno) ; 
if (yesno == 1 ) 
step = num_pts/cycles; 
else 
{ 
printf("\nEnter the desired Decimation. (e.g. for every point enter i,\n"); 
printf ("for every other point enter 2, for every third point enter 3, etc"); 
printf ("\nDecimation must be evenly divisible into the pts. per digit."); 
printf("\n(NOTE: Decimation must be less than *d tc ensure \n",num_pts); 
printf (" at least one point per digit is used!) \n"); 
printf ("\n Decimation: "); 
scanf ("%d", &step) ; 
if ( num_pts%step != 0) 
{ 
printf£("\nInvalid decimation. Try again."); 


printf ("\n Decimation: "); 
scanf ("%d", &step) ; 

. if ( num_ptststep != 0) 
print£("\nError. Bye!\n"); 
} 


if (step > num_pts) 
* { 
printf ("\nDecimation chosen is too large, use a smaller number.\n"); 











printf ("\n Decimation: "); 
scanf ("%d", &step) ; 
Lf (step > num_pts) 
print£("\nSorry! Still won't work, aborting\n"); 
} 
} 


/* Determine if coherent averaging of sequences is desired for improving 
processing gain. */ 


printf ("\nEnter the number of Sequence Frames to be coherently averaged."); 
printf ("\nNumber Frames = 1 implies no coherent averaging desired.\n\n"); 
printf (" Number Frames: "); 
scanf ("£c", 6num_coh_avg) ; 
if (num_coh_avg <= 0) 
{ 
printf ("\nInvalid Number of Frames. Must Use Positive Integer."); 
printf ("\n\n Number Frames: "); 
scanf ("&d", &num_coh_avg)? 
if (num_coh_avg <= 0) 
{ 
printf("\nInvalid Number of Frames. Aborting!\n"); 
exit (1); 
} 
} 


/* Query for Doppler range to be searched. An input of zero will cause 
the resampling for doppler steps to be skipped. Doppler is coverted 
from knots to meters/sec for computation. Sequence frequency resolution 
is computed to assist in determining doppler bins to be searched. */ 


printf("\nEnter the doppler range to be searched (knots)."); 
printf ("\n\n Doppler (+/-): "); 
scanf ("%£", doppler) ; 
printf ("\nOne knot corresponds to a %1.5f Hz shift", *fc-*fc*(1.0-0.5/snd_vel)); 
printf(" in frequnecy.\nEnter the Doppler increment to be searched."); 
printf ("\n\n Step (knots): "); 
scanf ("%f", &doppler_step) ; 
if (doppler_step == 0.0) 
doppler_step = 1.0; 
doppler /= 2.0; /* convert knots to meters/second */ 
doppler_step /= 2.0; 
doppler_start = -doppler; 
Goppler_end = doppler; 


/* Query for input and output files. */ 


print£("\nEnter Input File Name (20 Characters Maximum) .\n"); 
printf ("\n File Name: "); 
scanf ("ts", infile); 
printf ("\nEnter Output File Name (20 Characters Maximum) .\n"); 
printf ("\n File Name: "); 
scanf ("%s", outfile) ; 
if ((fpil = fopen(infile, "r") )==NULL) 
{ 
printf ("\nUnable to Open Input File. Aborting!\n"); 
exit (1); 
} 
if ((fp2 = fopen(outfile, "w") ) ==NULL) 
{ 
printf ("\nUnable to Open Output File. Aborting!\n"); 
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exit (1); 
) 


/* After Decimation there are less data pts => comput new num_pts. */ 
num_pts /= step; 


/* Set data length for input to maximum needed for highest doppler, 
which corresponds to the maximum number of points used in interpolation.*/ 


Ti = Ts/(1 + doppler_end/snd_vel); 
vect_len_ inter = Ti*vect_len/Ts + 1; 


/* Allocate memory as required. */ 


indatal = new(vect_len_ inter); 
indata2 = new(vect_len_inter); 
im_data = new(num_pts*seq_ len); 
re_data = new(num_pts*seq_len); 
mag = new(num_pts*seq len); 
phase = new(num_pts*seq_ len); 
avg_re = new(num_pts*seq_len); 
avg_im = new(num_pts*seq_len); 
h_data_im = new(seq_len+1); 
h_data_re = new(seq_len+1); 
data_i_re = new(vect_len); 
data_i_im = new(vect_len); 


/* DATA PROCESSING SECTION. 
Doppler processing is performed automatically according to the step 
size specified by the user. Data is processed in lengths corresponding 
one segment at a time. Sufficient points for complete interpolation 
of one doppler shifted sequence are picked off the input file. 
During doppler shift interpolation the input end point is saved for 
processing in the next segment to insure segment continuity. */ 


for (doppler = doppler start;doppler <= doppler_end; doppler+=doppler_step) 
{ 





/* Initialize demodulation routine. */ 
demodulate (fc, fs, doppler, vect_len) ; 


/* Determine the sampling period for the new sampling frequency used in 
interpolation, as a function of doppler. */ 
Ti = Ts/(1 + doppler/snd_vel); 


/* If zero doppler case is being performed endpoint bookkeeping is not 
used and data length is the same as input data (vect_len). */ 


if (doppler != 0 ) 
{ 
fscanf (fp1,"%f\n", Gindatal(0}); 
vect_len_inter = Ti*vect_len/Ts; 
indatal++; 
*indata2++ = 0; 
} 
else 
vect_len_inter = vect_len; 


/* Initialize flags and interpolation parameters. frirst_run indicates 








/* 


/* 


/* 





first segment being processed, and end indicates EOF reached. Time_s 
indicates current relative sampling time and time_i indicates current 
relative interpolation time. */ 


‘time_s = 0.0; 


time_i = Ti; 

end = FALSE; 

first_run = TRUE; a 
doppler_bin = PI * doppler/(2*snd_vel); 


Ta is sequence length in seconds. Header is printed to separate 
doppler scans. */ 


Ta = cycles * seq _len /(*fc * (1 + doppler/snd_vel)); 

fprintf (fp2,"START\n") ; 

fprintf(fp2,""'%6.2£ KNOTS DOPPLER'\n",2*doppler) ; 
fprintf(fp2,""SEQUENCE PERIOD (%7.4f SEC.) '\n",Ta); 

fprintf (fp2,"%5d %10.7f£\n",seq_len*num_pts, Ta/ (seq_len*num_pts-1)); 


Process entire input file for current doppler search. */ 


while (!end) 
{ 
for (j=0; j<vect_len_inter; j++) 
LE ((c#fgetc(fpl1) ) ==EOP) 
{ 
rewind(fpl) ; 


demodulate(fc,fs,dopplier,0); /* clears demodulator for next run. 
end = TRUE; 
break; 
} 
else 


{ 

ungetc(c, fpl); 
fscanf(fp1, "%f\n", Gindatal[j)); 
} 


*/ 


Test for EOF reached before complete segment read. Prevents processing 


partial segments. */ 


if (tend) 
{ 


Test for first segment. If it is the first data point has not been 
filtered, adjusts bookkeeping of endpoint. */ 


Lf ((doppler!=0.0) && first_run) 
{ 
--indatal; 
~~indata2; 
vect_len_inter += 1; 
} 
Lowpassing and then highpassing reduces to a BP equivelant and avoids 


complex numbers at this point. Caution: the filter phase must be linear 


in the pass region. */ 


hipass (indatal, indatai, vect_len_inter) ; 
lowpass3(indatal, indatal, vect_len_inter); 

demodulate (indatal, indata2,doppler_bin, vect_len_inter); 
lowpassi (indatal, indatal, vect_len_inter); 

lowpass2 (indata2, indata2,vect_len_inter); 
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/* This loop performs a linear interpolation of the demodulates. Maintains 
bookkeeping of endpoint, and first segment. Interpolates vect_len points 
for hadamard processing. Skips this if zero doppler case. */ 


if (doppler != 0.0) 
{ 
4=0; 
if (first_run) 
a { \ 
vect_len_inter -= 1; 
first_run = FALSE; 
) 
else 
{ 
-~-indatal; 
--indata2; 





} 
for (j=0; j<vect_len; j++) 


{ 


/* First case is down shift of doppler. xf] 


if (time_i > time_s+Ts) 
{ 
time_s += Ts; 
itt; 
h = time_i - time_s; 
data_i_im{j) = indatal(i) + (h/Ts)*(indatal (i+1)-indatal(i}); 
data_i_re{j} = indata2{i} + (h/Ts)*(indata2(it+1)-indata2[i}); 


/* Reset sampling time and interpolation time. The relative time is important 


. only. */ 
time_s = 0.0; 
time_i = h + Ti; 
Lf (time_i > 2*Ts) 
{ 
time_s += Ts; 
itt; 


} 
} 


/* Second case is upshift in doppler. */ 


else 

{ 
h = time_i - time_s; 
data_i_im(j] = indatai(i) + (h/Ts)*(indatal (it+1}-indatal({i}); 
data_i re(j) = indata2(i) + (h/Ts)*(indata2(itl]-indata2(i}); 
time_i += Ti; 

Lf (time_i > time_s + Ts) 
{ 
time_s += Ts; 
itt; 
time_i -= time_s; 
time_s = 0.0; 
} 








/* Save indatal and indata2 working space. */ 


temp _ptrl = indatal; 
temp_ptr2 = indata2; 
*temp_ptrl++ = indatai (iti); 
xtemp_ptr2++ = indata2{itl); 


/* Set pointers to interpolated values */ 
indatal = data_i_im; 


indata2 = data_i_re; 
} 


/* Decimate in time according to step. And compute sequence de level. */ 


n=0; 

dclvl_re = 0.0; 

delvi_im = 0.0; 

for (j=0;j<vect_len;j+=step) 
{ 
im_data[(n} = indatal(j)}; 
re_data(n++] = indata2(j); 
} 


/* Restore indatal and indata2 working space, when processing for doppler. 


if (doppler != 0.0) 
{ 
indatal = temp_ptrl; 
indata2 = temp_ptr2; 
} 


/* Scramble data vector and process data according to interleave via FHT 
(HADAMARD) and unscramble. Real and imaginary parts are processed 


separately. */ 


for (j=0; j<num_pts; j++) 
{ 
h_data_im[0} = 0; 
h_data_re[0} = 0; 
n=0; 


for (k=0; k<(seq_len*num_pts) ;k+=num_pts) 


{ 


h_data_im(scram(n]} = im_data(k+j); 


h_data_re(scram(n++]} = re_data[(k+j); 
} 

hadamard (degree, h_data_im); 

hadamard (degree, h_data_re) ; 

n=0; 


dclvl_re+=h_data_re[0); /* Save pedestal information. */ 


dclivl_im+=h_data_im{0); 


for (k=0;k<(seq_len*num_pts) 7k+=num_pts) 


{ 
im_data[k+j) = h_data_im{unscram[{n]); 


} 
) 








re_data[k+j]) = h_data_re[unscram(n++]); 


/* Compute pedestal corrections for real and imaginary parts. 


pdstal_re = (dclvl_re/num_pts)*rfact - (dclvl_im/num_pts) *ifact; 














*/ 





pdstal_im = (dclvl_re/num_pts)*ifact + (delvl_im/num_pts) *rfact; 


/* Case 1 is if Coherent averaging is performed. Else, no coherent averaging 
is performed. Each segment is output. xf 


if (num_coh_avg > 1) 
{ 


/* Make DC level correction. */ 


for (k=0;k<seq_len*num_pts;k++) 
{ 
avg_re{k}+= (re_data[{k}]-pdstal_re); 
avg_im{k)+= (im_data[k}-pdstal_im) ; 
} 
count++; 
if (count==num_coh_avg) 
{ 
count = 0; 
for (k=0;k<seq_len*num_pts; k++) 
{ 
avg_re(k)/=num_coh_avg; 
avg_im{k}/=num_coh_avg; 
} 


/* Compute magnitude and phase and print results. */ 


mag_phase(avg_re,avg_im,mag, phase, (num_pts*seq _len)); 
for (j=0; j<num_pts*seq_len; j++) 
{ 

fprintf(fp2,"%8.1f %8.1f\n",mag{j},phase{3j)); 
avg_re({j) = 0.0; 
avg_im({j) = 0.0; 
} 


/* Case 2. No coherent averaging. Compute magnitude and phase and 
print results. x/ 


else 
{ 
for (k=0;k<seqg_len*num_pts; k++) 
{ 
re_data(k]-= pdstal re; 
im_data(k)~-= pdstal_im; 
} 
mag_phase(re_data, im_data,mag, phase, (num_pts*seaq_len)); 
for (j=0; j<num_pts*seq_len; j++) 
fprintf£ (fp2,"%8.1f %8.1£\n",mag{j),phase{j)); 
} 
} 
} 
} 


/* Finished close out variables and files. */ 


fclose(fpi); 
f£close(fp2); 
free (fc); 
free (fs); 











free (indata}) ; 
free (indata2) ; 
free (im_data); 
free (re_data); 
free (mag); 
free (phase) ; - 
free (avg_re); 
free (avg_im) ; 
i free (h_data_im); 
i free (h_data_re); : 
free(data_i_re); 
free (data_i_im); 


} 
/* NEW is a short function to allocate memory for floating point vectors. */ 
float *new(size) 
int size; | 
{ 
float *newdata; 
if (( newdata = (float *)malloc(size*sizeof (float) ))==NULL) 
{ \ 
printf£("Cannot Allocate Storage!\n"); 
exit (1); | 
} 
| 


return (newdata) ; 
} 


B. INITIALIZATION PROGRAMS 
1. INIT_HAD 


This program prompts for the primitive polynomial initial register load and 
allocates memory space to the variables scram and unscram. Returns initial register 


load, initial. 

/* INIT_HAD is a program to initialize memory allocation for the scrambling 
and unscrambling arrays to be used with Past Hadamard Transform. It 
also returns the initial register load used with the Shift Register 
Generator. 

Variables: 
scram - external array to hold scrambling indices. 
unscram - external array to hold unscrambling indices. 
initial ~ initial register load for SRG. 


law - polynomial law that defines the SRG structure. 


length - length of sequence generated by the input law. Assumes 
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maximal length polynomial. 


temp - temporary holding spot for determining sequence length.*/ 


#include "macrofile.h" 
unsigned. init_had() 

{ 

extern unsigned law; 

unsigned initial, temp, length; 


printf ("\nEnter polynomial law for the desired M-Sequence. Use"); 
printf ("\noctal integer representation only."); 
printf("\n(e.g. D3 + D + 1 = 1011 binary => 13 Octal.)\n"); 
printf ("\n LAW: "); 
scanf ("%o", Glaw) ; 
temp = law; 
length = 1; 
while (temp>>=1) 
length<<=1; 
length--; 
Lf ((seram = (unsigned *)malloc(length*sizeof (unsigned))) == NULL) 
{ 
print£("\nCannot Allocate Scram Array!!!!\n"); 
exit (1); 
} 
if ({unscram = (unsigned *)malloc(length*sizeof(unsigned))) == NULL) 


{ 
printf ("\nCannot Allocate Unscram Array!!!!\n"); 





exit (1); 

} 
printf£("\nEnter the initial register load in Decimal, do not use zero."); 
printf ("\n\n Initial Lead: "); 


scanf ("$u", &initial) ; 
if (initial == 0) 
{ 
printf ("\nError! Initial load cannot be 0!!\n\n"); 
printf("\nEnter the initial register load in Decimal, do not use zero."); 
printf ("\n\n Initial Load: "); 
scanf ("%u", &initial); 
if (initial == 0) 
{ 
printf("\nError! Aborting.\n\n"); 
exit (1); 
} 
} 
return (initial); 


) 


2. FWD_HAD 
Computes the permutation indices for use with the FHT. Data is transmitted in 


the forward direction. Returns the degree of the primitive polynomial. 


/* FWD_HAD() computes the permutation indices for scrambling and 
unscrambling a data vector generated using Maximal Length Sequences. 
These indices are used with the Fast Hadamard Transform and is performed 
in place. Sequences are assumed to be transmitted in the foward 








direction. 
Variables: 
law - Shift Register Generator law to be used. 
initial - initial register law. 
scram - array of indices for scrambling data. 
unscram - restoration indices for use after FHT. 
degree - gives the degree of the polynomial law. 
temp ~ temporary holder for computing degree and seq_len. 
S_law - law formed from law to implement the S_gen structure. 


rev_law - law in the reverse order to implement the L_gen structure. 





end_bit - maintains end bit for around end feed in S gen. 
L_gen ~ variable that acts as L generator delay registers. 
S_gen - variable that acts as S generator delay register. 


Reference: The Feedback generators used for forming the foward scrambling 
and unscrambling indices are adapted from: 


M. Cohn and A Lempel, 'On Fast M-Sequence Transforms', 
IEEE Transactions on Information Theory, Jan. 1977. */ 


unsigned fwd_had(law, initial, scram,unscram) 
unsigned law, initial, *scram, *unscram; 

{ 

unsigned degree, temp, temp2,seq_len,S_ law, xrev_law; - 
unsigned end_bit,L_gen,S_gen; 5 


int i,j,count; 


/* Initialize variables. */ 
temp = law; 
degree = 0; 


S_gen = 0; 
rev_law = 0; 
seq_len = 1; 


/* Computes the length of the M~Sequence, and Polynomial Degree. */ 
while (temp>>=1) 

{ 

seg_len <<= 1; 

degreett; 


} 


/* Reverses law for use with L_gen. */ 
temp = law; 
for (1=0; i<= degree;it++) £ 
{ 
rev_law = (tempé1) | (rev_law<<1); 
temp>>=1; 
} a 














/* Set end_bit for register end feed for S_ gen logic. Set initial load 
for S_gen. Set law for use in S_gen logic.*/ 
end_bit = seq_len--; 
S_law = (law>>1); 
end_bit>>=1; 
’ temp = initial; 
for (i#0; i<seq_len-1; i++) 
{ 
if (tempé1) 
, temp = (temp*law)>>1; 
else 
temp >>= 1; 
if(i >= seq_len-degree) 
{ 
S_gen = S_gen} (tempé1); 
S_gen <<=1; 


} 
“S_gen J= (initialsé1); 
/* Load L_gen to generate unscrambling indices. */ 
L_gen = (1<<(degree-1)); 


/* Compute permutations scram and unscram. */ 
for (1=0; i<seq_len; i++) 
{ 
unscram(i} = L_gen; 
temp2 = 0; 
temp = S_gen; 
for (j=0; j<degree; j++) 
{ 
temp2 <<=1; 
temp2 |= (tempé&1); 
* temp >>=1; 
} 
scram{i] = temp2; 
if (L_gené1) 
7 L_gen = (L_gen*rev_law)>>1; 
else 
L_gen>>=1; 
temp = (S_gen & S law); 





/* Count the number of 1's for modulus two sum for end feedback to S_gen. */ 


count=0; 
for (j=0; j<degree; j++) 
{ 
if (tempé&1) 
{ 
countt+; 
temp>>=1; 
} 
else 
temp>>=1; 
} 
Lf ( (count%2) == 0) 
a S_gen <<= 1; 
else 
S_gen = (S_gen<<i) 11; 
S_gen &= seq_len; 
‘ } 
return (degree) ; 











3. REV_HAD 
Computes the permutation indices for use with the FHT. Data is transmittedin . 


the reverse direction. Returns the degree.of the primitive polynomial [Ref. 15]. 


/* REV_HAD() computes the scrambling and unscrambling indices for a 
Maximal Length Sequence that is to be processed using a Fast Hadamard 
Transform and is performed in place. The sequence is assumed to be 
transmitted in the reverse direction. 


law ~ The polynomial law defining the Maximal Length Sequence that is 
used. 


initial - The initial register load for the shift register generator. 


scram —- Pointer to the array to contain the scrambling indices, 
to be used with the transformed array. The transform 
is not done in place. 


unsccam ~ Pointer to the array to contain the unscrambling indices. 
The transform is not done in place. 


degree - The degree of the law is returned. 


Reference: This procedure was adapted from the article "On Fast 
M~Sequence Transforms", by Martin Cohn and Abraham Lempel, 
IEEE Transactions on Information Theory, Jan. 1977. 


The program was modified from code provided by Kirk 
Metzger with permission. xf 


finclude “macrofile.h" 


unsigned rev_had(law, initial, scram, unscram) 

unsigned int law, initial, *scram, *unscram; 

{ 

/* temp Scratch variable for use in finding degree etc. 
end bit Bit is set corresponding to highest bit in law. (i.e. the 

end register n). 

rev_initial Initial register load for reverse generation. 
seg_len Length of sequence based on law. 
index Loop counter. 
ss_contents Shift register contents for scram array. 
ms_contents M sequence register contents for unscram array.*/ 


unsigned int degree, temp, end_bit, rev_initial, seq_len, index; 
register unsigned ss_contents, ms_contents; 


/* Initialize variables. */ 
temp = law; 
seg_len = 1; 
rev_initial = 0; . 
degree = 0; 


/* Find sequence length and degree of polynomial. */ 
while (temp>>=#1) 
{ 
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seq_len<<=1; 
degreet+; 
} 


/* set end_bit */ 
end_bit = seq len-->>i; 


/* find reverse generator initial load */ 
for(index = 0; index < degree; index++) 
{ 
rev_initial = (rev_initial<<1) | (initialé1); 
initial>>=1; 
} 


/* generate scram and unscram values using law, end_bit and values of 
ss_contents and ms_contents. */ 
ms_contents = 1; 
ss_contents = rev_initial; 
for (index =0; index<seq_len; index++) 
{ 
*scramt+ = ss_ contents; 
*unscram++ = ms_contents; 
temp = ss_contentsélaw; 
ss_contents>>=1; 
do { if (temps1) 
ss_contents*send_bit; } 
while (temp>>=1) ; 
if (ms_contentsé1) 
ms_contents = (ms_contents“law) >>1; 
else 
ms_contents>>=1; 
} 
return (degree) ; 
} 


4. GtT_FLT_CO. 
Retrieves the filter coefficients to be used used in initializing the filtering 


programs. Two tenth order and one fifth order filters are used in this program. 


/* GET_FLT_COEF retrieves the filter coefficients from the wser specified 
file using a columnar format. The first column is the b coefficients 
and the second is the a coefficients. The first 11 rows are the highpass 
coefficients used for the upper end of the passband filter and the next 
eleven are the lowpass coefficients for the lower end of the passband. 
the final six rows are for the low pass filter used after demodulation. 


The coefficients file may be created by any program desirea by the user 
but must be in the appropriate format. However, a MATLAB program, 
set_filter.m, does this automatically. Set_filter.m uses a BUTTERWORTH 
filter (10th order) for the passband and a CHEBYCHEV filter for the 
lowpass filter (Sth order). */ 


finclude <stdio.h> 
finclude <math.h> 














void get_fit_coef (hi_a,hi_b, low_a, low_b, lowl_a,lowi_b) 
float *hi_a, *hi_b, *low_a, *low_b, *lLowl_a, *lowl_b; 

{ 

int i; 

FILE *fp; 

char filterfile{21)}; 


printf ("\nEnter the filename with the desired filter coefficients."); 
print£("(20 Characters Max) \n"); 
printf ("\n Filename: "); 
scarf ("%s", filterfile); 
if ((fp = fopen(filterfile, "r")) == -NULL) 
{ 
print£("\nUnable to open file. Aborting!\n"); 
exit (1); 
} 
for (i=0;i<=10;i++) 
fscanf(fp,"%f %£\n",&hi_b(i},éhi_aliJ); 
for (1=0;i<=10; i++) 
fscanf(fp,"%f %f£\n",&lowl_b(i],&lowl_al(i}); 
for (i=0; i<=5; 1++) 
fscanf(fp,"%f tf\n",&low_b(i)},&low_a[i)); 
fclose(fp); 
) 


€. DEMODULATION AND FILTERING 


1. DEMODULATE 
This function performs demodulation of the input carrier signal as well as the 
necessary frequency shift when compensating for Doppler. No assumption. are made 
about the sampling frequency within this program. It may be set and cleared as desired. 


The time index n is remembered so the repetitive calls may be made. 


/* DEMODULATE performs complex demodulation of a sinusoidal carrier for 
digital signal processing. Given an input signal, sampling frequency, 
and carrier freauency two outputs are produced, one tor the real part 
and one for the imaginary part. Tne user determines the number of points 
to be processed in each call to the program, 


sin_out - Input signal and output as imaginary part of the demodulated 
signal. Input data is overwritten and lost forever. 
Oa initialization sin_out is set to the desired carrier 
frequency. (floating type pointer.) 


cos_out - Output signal as real part of the demodulated signal. 
On Initialization cos_out is set to the desired sample frequency. 
(floating type pointer.) 


dopp bin - Denotes the doppler frequency shift in the sampling frequency 


to be used when doppler processing is used, and is determined 
externally. When no doppler processing is to be used it must 
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be set to 0. 
n - counter that maintains time increment for digit processing (static). 
theta - digital frequency (static). 


flag - Set when system has been initialized (static). Cleared when 
N has been set to zero. 


N - the number of data points to be processed. Whenever set to zero 
the system must be reset. 


Reset: call with all arguments, N=0. 
Set :First call with sin_out=fc, cos_out=fs, N is don't care. 


Operation : Call with sin_out as input data, returns sin_out as imaginary 
part and cos_out as real part. N= number of data points to 
process and is the same length as the input data. x/ 


#include <math.h> 
fdefine PI 3.1415926536 


void demodulate(sin_out,cos_out,dopp_bin,N) 
float *sin_out, *cos_out, dopp_bin; 

int N; 

{ 

int i; 

static int flag, n; 

static double theta; 


if (flag == 0) 
. { 
flag = 1; 
theta = 2 * PI *(*sin_out/(*cos_ out)); 
n= 0; 
return; 
} 
else if( N == 0) 
{ 
flag = 0; 
theta = 0.0; 
n= 0; 
return; 
} 
else 
{ 
for(i=0; i<N; i++) 
{ 
*cos_out++ = *sin out * cos(n*(theta + dopp bin)); 
xsin_out#+ = *sin_out * sin(n*(theta + dopp_ bin)); 
nt+4; 


} 
return; 
} 
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2. HIGHPASS 


This program performs a highpass filtering of input data. It is assumed to be 





used in conjunction with a lowpass filter to form a bandpass filter. Processing in this 
manner simplifies the processing by avoiding complex operations. It maintains its state 
on subsequent calls and can be initialized with any type filter up to order 20. It may 
also be cleared and reset. (Note: All filtering programs are generic and may be adapted 
for any filtering application, lowpass, highpass, bandpass and bandstop. The names | 


were chosen to distinguish them from each other within SEQREM.) 


/* HIGHPASS is a filtering program that operates using difference 
equations. It is intended that the filter be IIR but this is not 
required, (i.e. H(z) = B(z)/A(z) is a polynomial where the Bn's and 
An's specify the coefficients of the difference equation, which must 
of the same length. If they are not of the same order zeros must be 
appended.) 





Initialization: The first call sets the desired filter coefficients 
to be used throughout. No filtering is performed at this stage. 

In this case in_num specifies numerator and out_den specifies the 
denominator coefficients. M is the number of coefficients, ordered 
highest to lowest. (NOTE: out_den{0} is assumed to be one and is not 
actually set during the initialization.) 


Example 1: H(z) = 1/z; => M = 2, in_num(0) = 0, in_numf{1) = 1, 
out_den(0] = 1, out_den{ij = 0 


Example 2: H(z) = (0.52*2 + 1)/(z*2 + 0.5z + 0.6) => M = 3, 
in_num(0) = 0.5, in _num(1) = 0, in_num(2)} = 1, 
out_den[0) = 1, out_den(1) = 0.5, out_den(2) = 0.6 


Filtering: Subsequent calls perform the desired filtering. In this 
case in_num is the input data, out_den is the output data from the 
fiiter, M specifies the number of data points being processed, and 
may be of any length as long as sufficient space has been allocated 
to the input pointer and the output pointer. 


Clearing: To clear and reinitialize the filter after initializing once, 
simply set M to 0 and call again as explained above. 


Data Types: Inputs are pointers to array's of float, M is integer. 
Filter Order: Maximum order is 25, 


Variables: 


in_num ~ pointer to input numerater coefficients, or to 
input data segment. 











in_den - pointer to input denominator coefficients, or to 
output data segment. (NOTE: input and output segments 
may be the same, but input values will be overwritten.) 


A - Denominator filter coefficients (static). 

B - Numerator filter coefficients (static). 

Y - Output storage buffer used in recursion (static). 

X - Input storage buffer usec in recursion (static). 

N - Remembers filter order for computing output (static). 


M - Filter order during initialization, and number of points 
to be processed on subsequent calls. 


flag - maintains initialization status (static). 


Other Filters: LOWPASS1(), LOWPASS2(), and LOWPASS3() are coded 
identically to this program and work in the same 
manner. */ 


#include <math.h> 
#define MAXORD 25 
#define MAX 24 


void hipass (in_num, out_den,M) 

float *in_num, *out_den; 

int M; 

{ 

static float A(MAXORD}, B{MAXORD], Y{MAXORD], X{MAXORD]; 
static int flag, N; 

int i, 3; 


/* Loop clears output und coefficient values for reuse in new filter.*/ 
if (M==0) 
{ 
for (i=0;i<MAXORD; i++) 
{ 


A{i) = 0; 
Bli) = 0; 
Y{i) = 0; 
X{i) = 0; 
} 
flag = 0; 
return; 


} 
/* Loop sets coefficients of the desired filter on first entry if flag=0 */ 
/* Note that A[{0} is assumed 1, since it corresponds to the desired output*/ 
else if (flag == 0) 
{ 
B(0}) -* in_num[(0}; 
for (i=l; i<M;it+) 
{ 
B(i} = in_num(i); 
A(i) = out_den[i]; 
} 
flag = 1; 
N = M; 

















return; 
) 
/* Otherwise perform filtering operations. Output stored in out_den[). */ 
else 
{ 
for(j=0; j<M; j++) 
{ 
X(MAX} = in_num(3j); 
/* Start filteriing operation. */ 
Y(MAX] = B[O}] * X{MAX); 
for (i=1;i<N; i++) 
Y(MAX) = Y{MAX) + B(iJ}*X(MAX-i) - Afi} *Y(MAX-i); 
out_den[j) = Y{MAX); 
/* Perform time shift of data points stored in filter. */ 
for (i=0;i<MAX;it+) 
{ 
Y(i} = Y(i+1); 
X{i} = X(iti); 
} 
} 
return; 


) 


3. LOWPASS1 
This program performs a lowpass filtering of input data. It operates exactly 
the same as HIGHPASS. 


/* LOWPASS] is a filtering program that operates using difference 
equations. It is intended that the filter be IIR but this is not 
required. (i.e. H(z) = B(z)/A(z) is a polynomial where tne Bn's and 
An's specify the coefficients of the difference -quation, which must 
of the same length. If they are not of the same order zeros must be 
appended.) 


Initialization: The first call sets the desired filter coefficients 
to be used throughout. No filtering is performed at this stage. 

In this case in_num specifies numerator and out_den specifies the 
denominator coefficients. M is the number of coefficients, ordered 
highest to lowest. (NOTE: out_den[(0) is assumed to be one and is not 
actually set during the initialization.) 


Example 1: H(z) = 1/z; => M = 2, in_num(0) = 0, in_num(i) = 1, 
out_den{0} = 1, out_den[1) = 0 


Example 2: H(z) = (0.52°2 + 1)/(z*2 + 0.5z + 0.6) => M = 3, 
in_num(0) = 0.5, in_num{1} = 0, in_num[2) = 1, 
out_den[0] = 1, out_den{1} = 0.5, out_den{2] = 0.6 


Filtering: Subsequent calls perform the desired filtering. In this 
case in_num is the input data, out_den is the output data from the 
filter, M specifies the number of data points being processed, and 
may be of any length as long as sufficient space has been allocated 
to the input pointer ana the output pointer. 














Clearing: To clear and reinitialize the filter after initializing once, 
simply set M-to 0 and call again as explained above. 


Data Types: Inputs are pointers to array's of float, M is integer. 
. Filter Order: Maximum order is 25. 
Variables: 


in_num ~ pointer to input numerater coefficients, or to 
input -data segment. 


in_den - pointer to input denominator coefficients, or to 
output data segment. (NOTE: input and output segments 
may be the same, but input values will be overwritten.) 

A - Denominator filter coefficients (static). 

B - Numerator filter coefficients (static). 

Y - Output storage buffer used in recursion (static). 

X - Input storage buffer used in recursion (static). 


N - Remembers filter order for computing output (static). 


M~ Filter order during initialization, and number of points 
to be processed on subsequent calls. 


flag - maintains initialization status (static). 
Other Filters: HIGHPASS(), LOWPASS2(), and LOWPASS3() are coded 


» identically to this program and work in the same 
manner. x/ 





#include <math.h> 
#define MAXORD 25 
#define MAX 24 


void lowpassi (in_num, out_den,M) 

float *in_num, *out_den; 

int M; 

{ 

static float A[MAXORD}], B[MAXORD], Y(MAXORD], X[MAXORD]; 
Static int flag, N; 

int i,j 


/* Loop clears output and coefficient values for reuse in new filter.*/ 
L£ (M==0) 
{ 
for (is0;i<MAXORD;i++r) 
{ 
A[ij 
B(i} 
Y (i) 
X(i) 
} 
flag = 0; 
return; 
} 
/* Loop sets coefficients of the desired filter on first entry if flag=0 

















: /* Note that A[0] is assumed 1, since it corresponds to the desired output*/ 
5 else if (flag == 0) 
. { | 
BO} = in_num[0)}; 
for (Ll; i<M;i++) 
{ o 
B(i) = in_num[i); 
A(i}) = ont_den{i}; 
} 
flag = 1; > 
N = M; 
return; 
} 
/* Otherwise perform filtering operations. Output stored in out_den({). */ 
else 
{ 
for(j=0; 3<M; j++) 
{ 
X{MAX) = in_num[j); 
/* Start filteriing operation. */ 
Y{MAX] = B{O} * X(MAX}; 
for (181; i<N;4i++) 
Y(MAX) = Y{MAX) + B(i)*X(MAX-i) - A[i) *Y{MAX-i); 
out_den[(j) = Y(MAX); 
/* Perform time shift of data points stored in filter. */ 
for (1=0;i<MAX; i++) 
{ 
Y(i) = y[(itl); 
X(i} = X[(itl); 
} 
} 
return; 
} ‘ 


4. LOWPASS2 
This program performs a lowpass filtering of input data. It operates exactly 


the same as HIGHPASS. 


/* LOWPASS2 is a filtering program that operates using difference 
equations. It is intended that the filter be IIR but this is not 
required. (i.e. H(z) = B(z)/A(z) is a polynomial where the Bn's and 
An's specify the coefficients of the difference equation, which must 
of the same length. If they are not of the same order zeros must be 
appended.) 


Initialization: The first call sets the desired filter coefficients 
to be used throughout. No filtering is performed at this stage. 

In this case in_num specifies numerator and out_den specifies the 
denominator coefficients. M is the number of coefficients, ordered 
highest to lowest. (NOTE: out_den[0} is assumed to be one and is not 
actually set during the initialization.) 





Example 1: H(z) = 1/z; => M = 2, in_num{0} = 0, in_num{1J = 1, 
out_den(0) = 1, out_den{1} = 0 











Example 2: H(z) = (0.5z*2 + 1)/(2*2 + 0.52 + 0.6) => M = 3, 
in_num{0] = 0.5, in_num{1) = 0, in_num[2) = 1, 
out_den(0] = 1, out_den({1] = 0.5, out_den[2] = 0.6 


> Filtering: Subsequent calls perform the desired filtering. In this 
case in_num is the input data, out_den is the output data from the 
filter, M specifies the number of data points being processed, and 
may be of any length as long as sufficient space has been allocated 

; to the input pointer and the output pointer. 


Clearing: To clear and reinitialize the filter after initializing once, 
simply set M to O and call again as explained above. 


Data Types: Inputs are pointers to array's of float, M is integer. 
Filter Order: Maximum order is 25. 
Variables: 


in_num - pointer to input numerater coefficients, or to 
input data segment. 


in_den - pointer to input denominator coefficients, or to 
output data segment. (NOTE: input and output segments 
may be the same, but input values will be overwritten.) 


A ~ Denominator filter coefficients (static). 


B - Numerator filter coefficients (static). 


«K 
1 


Output storage buffer used in recursion (static). 


X - Input storage buffer used in recursion (static). 


=z 
1 


Remembers filter order for computing output (static). 


M - “.lter order during initialization, and number of points 
to be processed on subsequent calls. 





flag - maintains initialization status (static). 


Other Filters: HIGHPASS(), LOWPASS1(), and LOWPASS3() are coded 
identically to this program and work in the same 
manner. xf 


#include <math.h> 
#define MAXORD 25 
#define MAX 24 


void lowpass2 (in_num, out_den,™) 

float *in_num, *out_den; 

int M; 

{ 

static float A[MAXORD], B[MAXORD}, Y(MAXORD], X[{MAXORD]; 
. static int flag, N; 

int i,j; 


/* Loop clears output and coefficient values for reuse in new filter.*/ 
. if (M==0) 
{ 























for (1=0;i<MAXORD; i++) 
{ 
A{i] = 0; 
B[i} = 0; 
Y{i] + 0; 
X{i] = 0; 
} 
flag = 0; y 
return; 
} 
/* Loop sets coefficients of the desired filter on first entry if flag=0 */ 
/* Note that A[0) is assumed 1, since it corresponds to the desired output*/ 
else if (flag == 0) 
{ 
B(O) = in_num[0}; 
for (i=1; i<M;itt) 
{ 
B(i) = in_num(i); 
Afi) = out_denf{i]; 
) 
flag = 1; 
N = M; 
return; 
} 
/* Otherwise perform filtering operations. Output stored in out_den{). */ 
else 
{ 
for(j=0; j<M; j++) 
{ 
X(MAX] = in_num(3); 
/* Start filteriing operation. */ 
Y{MAX} = B(O) * X{MAX]; 
for (i=1;1<N; i++) 
Y(MAX) = Y{(MAX) + B(L]*X(MAX-i) - A(i]*¥ (MAX~i); 
out_den[{j] = Y(MAX]; 
/* Perform time shift of data points stored in filter. */ 
for (1=0;i<MAX; i++) 
{ 
YCi) = Y(i+1); 
X(i} = X{it1j; 
} 
} 
return; 


) 


5. LOWPASS3 
This program performs a lowpass filtering of input data. It operates exactly 
the same as HIGHPASS. 


/* LOWPASS3 is a filtering program that operates using difference 
equations. It is intended that the filter be IIR but this is not 
required. (i.e. H(z) = B(z)/A(z) is a polynomial where the Bn's and 
An's specizy the coefficients of the difference equation, which must 
of the same length. If they are not of the same order zeros must be 
appended.) 
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Initialization: The first call sets the desired filter coefficients 

to be used throughout. No filtering is performed at this stage. 

In this case in_num specifies numerator and out_den specifies the 

denominator coefficients. M is the number of coefficients, ordered | 
y) highest to lowest. (NOTE: out_den[{0}) is assumed to be one and is not 

actually set during the initialization.) 


Example i: H(z) = 1/z; => M = 2, in_num[0) = 0, in_num{i) = 1, 
, out_den(0} = 1, out_den{1) = 0 


Example 2: H(z) = (0.52*2 + 1)/(z*2 + 0.5z + 0.6) => M = 3, 


in_num{0) = 0.5, in_num[(1} = 0, in_num(2) = 1, 
out_den[0) = 1, out_den{1) = 0.5, out_den[2] = 0.6 


Filtering: Subsequent calls perform the desired filtering. In this 
case in_num is the input data, out_den is the output data from the 
filter, M specifies the number of data points being processed, and 
may be of any length as long as sufficient space has been allocated 
to the input pointer and the output pointer. 


Clearing: To clear and reinitialize the filter after initializing once, 
simply set M to 0 and call again as explained above. 


Data Types: Inputs are pointers to array's of float, M is integer. 
Filter Order: Maximum order is 25. 
Variables: 


in_num - pointer to input numerater coefficients, or to 
. input data segment. 


in_den - pointer to input denominator coefficients, or to 
output data segment. (NOTE: input and output segments 
may be the same, but input values will be overwritten.) 

A ~ Denominator filter coefficients (static). 

B - Numerator filter coefficients (static). 

Y - Output storage buffer used in recursion (static). 

X - Input storage buffer used in recursion (static). 


N - Remembers filter order for computing output (static). 


M - Filter order during initialization, and number of points 
to be processed on subsequent calls. 


flag - maintains initialization status (static). 
Other Filters: HIGHPASS(), LOWPASS1(), and LOWPASS2() are coded 


identically to this program and work in the same 
‘ manner. */ 


#include <math.h> 
define MAXORD 25 
. #define MAX 24 











void lowpass3(in_num,out_den,M) 

float *in_num, *out_den; 

int M; 

{ 

static float A([MAXORD]), BIMAXORD], Y(MAXORD], X{MAXORD); 
static int flag, N; 

int i,47 


/* Loop clears output and coefficient values for reuse in new filter.*/ 
if (M==0) 
{ 
for (i=0;i<MAXORD; i++) 
{ 
A({i] = 0; 
B{i} 
Y{i} = 0; 
X(i] = 0; 
} 
flag = 0; 
return; 
} 
/* Loop sets coefficients of the desired filter on first entry if flag=0 */ 
/* Note that A({0} is assumed 1, since it corresponds to the desired output*/ 
else if (flag == 0) 
{ 
B(O) = in_num[(0}; 
for (i=l; i<M;it++) 
{ 
B{ij} = in_num(i}; 
A(i} = ovt_den[i]; 
} 
flag = 1; 
N= M; 
return; 
} 
/* Otherwise perform filtering operations. Output stored in out_den[]. */ 
else 
{ 
for(j=0; j<M; j++) 
{ 
X(MAX} = in_num(J]; 
/* Start filteriing operation. */ 
Y(MAX) = B[O) * X{MAX]; 
for (i=1;i<N;i++) 
Y(MAX] = Y{MAX]) + BCL) *X(MAX-i} - Afi) *Y(MAX-ij; 
out_den(j) = Y(MAX); 
/* Perform time shift of data points stored in filter. */ 
for (i=0;1i<MAX; i++) 
{ 
Yfi) = Y(itl); 
Xfi) = X(itl); 
} 
} 
return; 


} 
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D. FAST HADAMARD TRANSFORM (FHT) 
1. HADAMARD 
This function performs the FHT of an input data vector in place. A power of 
two length is required. Computation is identical to the FFT butterfly process with the 


complex mul-iplies set to one. 


/* HADAMARD performs the Fast Hadamard Transform (FHT) on the input data. 
The transform is performed in place and assumes no interleaving. 
Interleaving must be performed externally. 


The FHT performs the same operations as the FFT with the complex 
multiplies set to + or - 1. The input and output data vector 
must be first ordered and then reordered after exit. 


Variables: 


degree - Provides the degree of the polynomial being worked with 
which gives the length of the M-Sequence. 


din - Pointer to the input data to be transformed. On exit it holds 
the transformed data. 


m,n - counters. 


itop, ispace, iwidth,ibot - common place holders used performing 
the butterfly operations of the FFT. */ 


#include <math.h> 


void hadamard (decree, din) 

float *din; 

int degree; 

{ 

int m,n, itop,ispace, iwidth, ibot; 
float temp; 

unsigned int one, size; 


one = ]; 
size = (one<<degree) ; 
for ( m@=1; m<=degree;mt+) 
{ 
ispace = (one<<m); 
iwidth = (one<<(m-1)); 
for (n=0;n<iwidth;n++) 
for( itop=n; itop<=«(size-2) ;itopt=ispace) 
{ 
ibot = itop + iwidth; 
temp = dinfibot]; 
din(ibot) = din(itop} - temp; 
din{itop) = din[itop] + temp; 
} 
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E. MAGNITUDE AND PHASE 
1. MAG PHASE 
This function simply computes the magnitude and phase of the demodulates =< 
after either coherent averaging or directly after the FHT and the appropriate 


permutations. 


/* MAG PHASE computes the magnitude and phase of the input data, and 
returns the results in the pointers maq and phase. 


realpart - Pointer for the real part of the input data of length N. 
impart - Pointer for the imaginary part of the input data of length N. 


mag - Pointer to data segment to hold the magnitude of the data and 
is also of length N. 


phase - Pointer for data segment to hold the phase of the data and is 
also of length N. (Degrees) 


N - Integer indicating the size of the input data segment. */ 


Sinclude <math.h> 
&define PI 3.1415927 


void mag_phase (realpart, impart, mag, phase, N) 
float *realpart, *impart, *mag, *phase; 
int N; 

{ 


int i; 


fox (i=s0;i<N; i++) 

{ 

mag(i) = sqrt (realpart(i}) * realpart[{i] + impart(i) * impart(il); 

if (realpart{i! == 0.0) 
{ 
puase[i}] = 0.0; 
} 

else 
phase(i} = atan(impart{i} / realpart{i})) * (180.0/PI); 


} 
} 











F. UTILITY PROGRAMS | 
1. SET_FILTER 
This is a handy MATLAB pregram for finding the coefficients for a tenth | 

order Butterworth-bandpass filter and a fifth order Chebychey lowpass filter, which | 


stores them in a file formatted for use with GET_FLT_COEF. Sampling frequency 





and cutoff frequencies are prompted from the user. 


% function set_filter Set_filter computes the filter coefficients for 

% the bandpass filter (thigh and lowpass filters) 

% and for the lowpass filter used in the seguence 
% removal program. The Dandpass region uses a 

% 10th order BUTTERWORTH { lter and the lckpass region J 
% uses a 5th order CHEBYCHEV filter with the results 

% stored in the user specified file in the format 

% required by the GET_FLT_CQEF function. The user is 
% queried for cutoff frequencies and samoling 

% frequency. 


function set_filter 





fprint£('\n WARNING! ifil! \n'); 
fprintf£('\n File to hold results should be non-existant prior to running\n'i: 
fprintf(' this routine because the results are appendec to any vcre-existing\n'); 
fprint£(' data.\n'); 
filename = input('Enter filename to store results: ','s'); 
fs = input('Ester the Sampling frequency to be used: ");7 
fc = input('Enter the Cut Off Frequency for the Low Pass Filter: '); “| 
fprintfé('\n Enter the Cut Off Freauencies for the BP filter.\n'); 
fl = input('Lower: '); 
fh = input('Upper: '); 
fs = f£s/2; 
iz (f1 > fs) [{fe > fs) 1 (fh > £s) 

error (‘Invalid Frequencies. Cutoft Frequencies Must be <= Fs/2"'); 
end | 
%& Finds the Coefficients for the High Pass Filter (BP low end} 
wn = £1/fs; 
(b1,a1] = butrer(10,wn, "high'); 
% Finds the Coefficients for the Low Pass Filter (BP hi end) 
wn = fh/fs; H 
(b2,a2} = butter (10,wn); 
% Finds the Coefficients for the Low Pass Filter. 
e wn = fc/fs; 
{b3,2a3] = cheby1(5,0.5,wn); 
fox i=1:11 

fprintf(filename,'%te t%e\n',bi(i),al(i)): 
end 
for isl:11 














fprintf(filename,'%te %e\n',b2(i),a2(i)); 
end 
for 11:6 

fprintf(filename,‘%e %e\n',b3(1),a3(4)); 
end 





2. APLOT 
This program was used to generate the plots in Appendix A with NCAR 


. 


graphics. It is data specific because of the nature of NCAR graphics when used for 


surface plotting routines. 


Cc This is a program to plot the results of the program SEQ REM 
c in successive plots for each doppler searched. 


INTEGER X,Y 
PARAMETER (X=510, Y=30) 


CHARACTER*32 IFILE, TEMP, SEQLEN 

CHARACTER*20 DOPPLER 

INTEGER I,J, LENGTH, ISZ 

REAL MAG(X,Y),PH(X,Y) , XPOS, YPOS, ZPOS 

REAL WORK (2*X*Y+X+Y) 

REAL TS,ANGH, ANGV, ZMAX 

COMMON /SRFIP1/ IFR, ISTP, IROTS, IDRX, IDRY, IDRZ, IUPPER, ISKIRT, 
1 NCLA, THETA, HSKIRT, CHI, CLO, CINC, ISPVAL 


35 FORMAT (A) 
40 FORMAT (A20) 


IFR = 0 

IDRZ = 1 

WRITE (*,*) ‘THIS PLOTS THE RESULTS FROM SEQREM PROGRAM' 
WRITE (*,*)-* ENTER THE INPUT FILE: ' 
READ (*, 40) IFILE 

OPEN (33, FILE=IFILE, STATUS='OLD') 
ANGH = -60. 

ANGV = 30. 

ZMAX = 0. 

CALL GOPKS (6, ISZ2) 

CALL GOPWK(1,2,1) 

CALL GACWK (1) 

READ (33, 35) TEMP 


50 READ (UNIT=33, FMT=*) DOPPLER 
READ (UNIT=33, FMT=*) SEQLEN 
READ (UNIT=33, FMT=*) LENGTH, TS 
Jn] 


60 d= J+ 
DO 80 I=1, LENGTH 
READ (UNIT=33, FMT=*, ERR=100, END=200) MAG(I, J), PH(I,d) 
IF (MAG(I,J) .GI. 2MAX) THEN 
ZMAX = MAG (I,J) 
ENDIF 
80 CONTINUE 





a, 








100 


200 





GO TO 60 


J = 5-1 

CALL GSELNT(0) 

CALL GSTXAL (2, 3) 

CALL GSCHH(.02) 

CALL GTX(0.5,0.9, ‘ARRIVAL TIME PLOT ') 
CALL GSTXAL (2, 3) 

CALL GTX (0.5, 0.85, DOPPLER) 

CALL EZSRFC (MAG, LENGTH, J, ANGH, ANGV, WORK) 
XPOS = 0. 

YPOS = -0.25*ZMAX 

ZPOS = 0.0 

CALL PWRZS (XPOS, YPOS, ZPOS, SEQLEN, 31, 25, 1,3, 0) 
CALL FRAME 

ZMAX = 0, 

GO TO 50 


J = 3-1 

CALL GSELNT (0) 

CALL GSTXAL (2, 3) 

CALL GSCHH (.02) 

CALL GTX(0.5,0.9, ‘ARRIVAL TIME PLOT *) 
CALL GSTXAL (2, 3) 

CALL GTX (0.5, 0.85, DOPPLER) 

CALL EZSRFC(MAG, LENGTH, J, ANGH, ANGV, WORK) 
XPOS = 0. 

YPOS = -0.25*ZMAX 

2POS = 0.0 

CALL PWRZS (XPOS, YPOS, ZPOS, SEQLEN, 31, 25, 1,3, 0) 
CALL FRAME 

CLOSE (33) 

CALL GDAWK (1) 

CALL GCLWK (1) 

CALL GCLKS 

WRITE (*, *) "FINISHED! 

STOP 

END 
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3. DOPPLOT 


This program was used to generate the. arrival time plot, Figure 3.1, with 


NCAR graphics. It is data specific because of the nature of NCAR graphics when used 


for surface plotting routines. é 
Cc This is a program to plot the results of the program SEQ_REM 
c making a doppler versus time plot for each doppler search. 


INTEGER X,Y 
PARAMETER (X=255, Y=23) 


CHARACTER*32 IFILE, TEMP, SEQLEN 

CHARACTER*20 DOPPLER 

INTEGER I,J, LENGTH, ISZ 

REAL MAG(X,Y),PH(X,Y),XPOS, YPOS, Z20S, TEMP] (X, Y) , TEMP2 (X, Y) 
REAL WORK (2*X*Y+X+Y) 

REAL TS,ANGH, ANGV, ZMAX 

COMMON /SRFIP1/ IFR, ISTP, IROTS, IDRX, IDRY, IDRZ, IUPPER, ISKIRT, 
1 NCLA, THETA, HSKIRT, CHI, CLO, CINC, ISPVAL 


35 FORMAT (A) 
40 FORMAT (A20) 


IFR = 0 
IDRZ = 1 | 
WRITE (*,*) THIS PLOTS THE RESULTS FROM SEQREM PROGRAM IN DOPPLER! ‘ 
WRITE (*,*) "VS, TIME FORMAT! | 
WRITE(?,*) "ENTER THE INPUT FILE: ' | 
READ (*, 40) IFILE 
OPEN (33, FILESIFILE, STATUS='OLD') : 
ANGH = -60. | 
ANGV = 30. 

2MAX = 0. 

CALL GOPKS (6, ISZ) 

CALL GOPWK(1, 2,1) 

CALL GACWK (1) 

READ (33, 35) TEMP 

Jnl 





50 READ (UNIT=33, FMT=*) DOPPLER 
READ (ONIT=33, FMT=*) SEQLEN 
READ (UNIT=33, FMI=*) LENGTH, TS 
Keil 
DO 53 I=1, LENGTH 
READ (33, *) TEMP1 (I, K) , TEMP2 (I, K) 
53 CONTINUE 
60 J = J+i 
DO 80 I=1, LENGTH 
READ (UNIT=33, FMT=*, ERR=50, END=200) MAG(I,J) , PH(I,d) 
IF (MAG(I,J) .GT. ZMAX) THEN 
2MAX = MAG(I,J) 
ENDIF 
80 CONTINUE o 
82 CONTINUE 


ra 
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DO 85 I=], LENGTH 
READ (UNIT=33, FMT=*, ERR=50, END=200) TEMP1 (I,K) , TEMP2 (I, K) 
85 CONTINUE 
GO TO 82 


200 Js J-1 
CALL GSELNT (0) 
CALL GSTXAL (2, 3) 
CALL GSCHH(.02) 
CALL GTX(0.5,0.9, ‘ARRIVAL TIME VS DOPPLER‘) 
CALL GSTXAL(2, 3) 
CALL GTX(0.5,0.85,'=5 TO 5 KNOTS DOPPLER’) 
CALL EZSRFC (MAG, LENGTH, J, ANGH, ANGV, WORK) 
XPOS = 0. 
YPOS = -0.25*ZMAX 
2POS = 0.0 
CALL PWRZS (XPOS, YPOS, ZPOS, "ARRIVAL TIME 22.36 SEC.',31,25,1,3,0) 
CALL FRAME 
CLOSE (33) 
CALL GDAWK(1) 
CALL GCLWK (1) 
CALL GCLKS 
WRITE (*, *) 'PINISHED' 
STOP 
END 


4. MGEN 


MGEN generates a file with a m-sequence in the forward direction in column 
form, given a primitive polynomial. The sequence is in (1,-1) form. A second file is 
also created, which holds one period of the register states. This can be useful for 


debugging and-as an aid for understanding m-sequences. 


/* This is a generator routine to generate the output from the least 
significant register in a M-Sequence shifter register generator. 
Output is in the form of -1 and 1 where the mapping is {-1,1} -> 
{1,0}. Register states are also stored to a specified file in decimal 
form. 


Input: law = the polynomial law specifying the M~Sequence to be 
generated. 
initial = the initial load placed in the Shift Registers. 


Output: output is a column format of -1 and 1's in floating format 
form for ease of use. output 2 is a column of numbers in decimal 
representing the register states over one period. 


Initiation: use 'MGEN' 


x] 

#include <malloc.h> 
#include <stdio.h> 
#include <math.h> 











main () 

{ 

float *mout; 

unsigned int temp, initial, law, seq_len; 
int i; 

char filename(20); 

FILE *fpl,*fp2; 





printf("\nEnter the Polynomial Law in Octal Form: "); 
scanf ("%o", &law) ; 
printf ("\nEnter the Initial Register Load in Octal(Not Zero!): "); 
scanf ("to",éinitial) ; 
if (initial==0) 
{ 
printf£("\n0 is Invalid Register Load. Aborting!\n"); 
exit (1); 
} 
temp = law; 
seq_len = 1; 
while (temp>>=1) 
seq_len<<=1; 
seq_len--; 
if ((mout = (float *)malloc(seq_len*sizeof (mout) ))==NULL} 
{ 
print£("\nCannot Allocate Memory to Output Variable\n"); 
exit(1); 
} 
printf ("\nEnter file to store M-sequence. (20 Characters max."); 
scanf ("%s", filename) ; 
fpi= fopen(filename, "w"); 
printf("\nEnter fire to store register states. (20 Characters max."); 
scanf("%s", filename) ; 
fp2= fopen(filename, "w") ; 
for (i=0;i<seq_ len; i++) 
{ 
fprintf(fp2,"su\n", initial) ; 
if (initials) 
{ 
mout (ij = -1.0; 
initial = (initial*law) >>1; 
} 
else 
{ 
mout (ij = 1.0; 
initial>>=1; 
} 
} 
fclose (f£p2); 
for (i=0;i<seq_len; i++) 
fprintf£ (fpl,"%£\n", *moutt+) ; 
fclose(fpl); 
} 
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5. MAKEFILE 
This is a file which makes compiling and linking of C programs easier on 
UNIX based systems. Files in the list are compiled, if not already, and are linked 
automatically with the MAKE command. The executable code is stored in the file 


following the -o flag, and -Im includes the necessary library routines. 


#makefile for thesis work. 


OBJECTS = seqrem.o rev_had.o init_had.o fwd_had.o get_flt_coef.o low_filterl.o low_filter2.o 
low_filter3.o hi_filter.o demod.o hadamard.o mag_phase.o 


SOURCES = seqrem.c set_had.c init_had.c fwd_had.c get_flt_coef.c low_filterl,> low ‘ilter2.c 
low_filter3.c hi_filter.c demod.c hadamard.c mag_phase.c 


LIBS = -im 
CFLAGS = -O 
seqrem: $ (OBJECTS) 
cc $(CFLAGS) S$(OBJECTS) -o seqrem $(LIBS) 


$ (OBJECTS) : macrofile.h 


6. EXAMPLE FILTER COEFFICIENT FILE 
This shows the structure of the format used to store the filter coefficients used 
in the various filtering programs. The numerator coefficients constitute the first 
column and the denominator the second. The first eleven rows are for the first tenth 
order filter, the next the second, and finally the last six make up the fifth order lowpass 


filter. 
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rant Sars 





6.127443e-03 
-6.127443e-02 
2.757350e-01 
-7,352932e-03 
1.286763e+00 
-1.544116e+00 
1.286763e+00 
-7,352932e~01 
2.757350e-01 
-6.127443e-02 
6.127443e-03 
6.127443e-03 
6.127443e-02 
%.757350e-01 
7.352932e-01 
1.286763e+00 
1.544116e+00 
1.286763e+00 
7.352932e-01 
2..757350e-01 
6.127443e-02 
6.127443e-03 
.438164e-05 
.190822e-05 
-438164e~-04 
.438164e-04 
.190822e-05 
.438164e-08 


1 
7 
1 
i 
7 
1 


1.000000e+00 


-9.960099e-01 


1.759559e+00 


-1.112098e+00 


8.747430e-01 


-3.474470e-01 


1.441634e-01 


~3.307116e-02 


6.691890e-03 


~6.798897e-04 


3.914710e-05 
1.000000e+00 
9.960099e-01 
1.759559e+00 
1.112098e+00 
8.747430e-01 
3.474470e-01 
1.441634e~02 
3.307116e~-02 
6.691890e-03 
6.798897e-04 
3.924710e-05 
1.000000e+00 
~4.515326e+00 
8.278174e+00 
-7.695624e+00 
3.625292e+00 
-6.920557e-01 
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